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ON RELIABLE CONTROL SYSTEM DESIGNS 
by 

JOHN DOUGLAS BIKDWELL 
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for the Degree of Doctor of Philosophy. 

ABSTRACT 


This report contains a method of approach and theoretical 
framework which advances the state of the art in the design of 
reliable multivariable control systems, with special emphasis on 
actuator failures and necessary actuator redundancy levels. 

The mathematical model consists of a linear time invariant 
discrete time dynamical system. Configuration changes in the 
system dynamics, (such as actuator failures, repairs, introduction 
of a back up actuator) are governed by a Markov chain that includes 
transition probabilities from one configuration state to another. 
The performance index is a standard quadratic cost functional, 
over an infinite time interval. 

If the dynamic system contains either process white noise 
and/or noisy measurements of the state, then the stochastic 
optimal control problem reduces, in general, to a dual problem, 
and no analytical or efficient algorithmic solution is possible. 
Thus, the results are obtained under the assumption of full state 
variable measurements, and in the absence of additive process 
white noise.- 

Under the above assumptions, the optimal stochastic control 
solution can be obtained. The actual system configuration can 
be deduced with an one step delay. The calculation of. the optimal 
control law requires the solution of a set of highly coupled 
Riccati-like matrix difference equations; if these converge (as 
the terminal time goes to infinity) one has a reliable design with 
switching feedback gains, and, if they diverge, the design is 
unreliable and the system cannot be stabilized unless more reliable 
actuators or more redundant actuators are employed. For the 
reliable designs, the feedback system requires a switching gain 
solution, that is, whenever a system change is detected, the feed- 
back gains must be reconfigured. On the other hand, the necessary 
reconfiguration gains can be precomputed, from the off-li: 
tions of the Riccati-like matrix difference equations. 
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Through the use of the matrix discrete minimum principle, a 
suboptimal solution can also be obtained. In this approach, one 
wishes to avoid the reconfiguration of the feedback system, and 
one wishes to know whether or not it is possible to stabilize the 
system with a constant feedback gain, which does not change even ^ 
if the system changes. Once more this can be deduced from another" 
set of coupled Riccati-like matrix difference equations. If they 
diverge as the terminal time goes to infinity, then a constant 
gain implementation is unreliable, because it cannot stabilize' the' 
system. If, on the other hand, there exists an asymptotic solution 
to this set of Riccati-like equations then a reliable control 
system without feedback r e c onf igurat i on can be obtained. The 
implementation requires constant gain state variable feedback, an& 
the feedback gains can be calculated off-line. 

In summary, these results can be used for off-line studies- 
relating the open loop dynamics, required performance, actuator 
mean time to failure, and functional or identical actuator 
redundancy, with and without feedback gain reconfiguration 
strategies. 


Thesis Supervisor: Michael Athens 

Title: Professor of Electrical Engineering and Computer Science 
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INTRODUCTION 

1 . 1 Motivation for the Research . 

This report addresses some of the current problems in interfacing 
systems theory and reliability, and puts this research in perspective 
with the open questions in this field. Reliability is a relative concept; 
it is, roughly, the probability that a system will perform according 
to specifications for a given amount of time. The motivating question 
behind this report is: What constitutes a reliable system? 

Knowledge of the reliability of a system is crucial. In this 
report, a system is reliable if it has a (quantitative) reliability of 
one, i.e., if the probability that the system will not perform according 
to specifications for a given period of time is zero. Therefore, the 
question "What constitutes a reliable system?" can be restated as: 

What are the specifications which a system must meet in order to be 
reliable? 

A system is normally designed in two stages: First, the components 

are selected in such a way as to meet the relicibility specifications; 
second, the control problem is formulated and solved for that configure- . 
tion of components. Although this procedure is over-simplified, it 
illustrates a second question: Should the control problem influence the 

choice of the configuration, and if so, how can this be achieved? The 
first part of the question is answered by history: The control problem 

influences configuration design now by iteration between the two stages 
of design. This is most likely not the best method! If a theory were 
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available which allowed a comparison between alternate designs, based 
on both the expected system reliability and the expected system perfor- 
mance, it would greatly simplify the current design methodology. It is 
unfortunate that at present there is no accepted methodology for a 
determination of expected system performance which accounts for changes 
in the performance characteristics due to failure, repair or reconfigura- 
tion of system functions. This report presents such a methodology for a 
specific class of linear systems with quadratic cost criteria. 

1.2 General Nature of the Problem . 

This Section presents the general theoretical framework necessary to 
approach the problem of reliable control system design. First, a 
discussion of scmie of the concepts in reliability theory will be. present- 
ed. The control- theoretic framework for the specific topics covered in 
this report will then be developed. Finally, the interrelationships- 
between systems theory and reliability theory will be explored, leading 
to a mathematical formulation of the reliable control system design . 
problem and a discussion of the general nature of the results presented 
in the remainder of this report. 

1.2.1 Reliability Theory. 

The generally accepted definition of reliability is stated in 
Appendix 1. Basically, the reliability of a system is the probability 
that the system will perform according to specifications for a given 
amount of time. In a system- theoretic context, the specification which 
a system must meet is stability; also, since, at least for most mathemati- 


cal models of systems, stability is a long-term attribute of the system. 
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the amount of time for which the system must remain stable is taken -to 
be infinite. Therefore, the following definitions of system reliability 

’f 

are used in this report: 

Definitioni 1 : A system (implying the hardware configuration, or mathe- 

maticaiL 'model of that configuration, and .its associated control and 
estimation' structure) has reliability r where r is the probability that 
the system* .will be stable for all time. 

Definition 2 : A system is said to be reliable if r = 1. 

Definition 3 ; A system design, or configuration, is reliable if it 

C “ ■ 

is stabilizable with probability one. 

These- "definitions of reliability depend on the definition of stability, 
and for*! -systems which can have more than one mode of operation, stability 
is not-' that easy to determine. In this report, stability will mean 
either"* mean-square stability (over some random space which will be left 
unspecified for the moment) , or cost-stability (again, an expectation 
over 'a- certain random space ) , which is basically the property that the 
accumulated cost of system- operation is bounded with probability one. 

(The definition of cost is also deferred.) 

The reliability of a system will depend on the reliabilities of its 
various components and on their interconnections. Thus, the systems 
engineer must have an understanding of the probabilistic mechanisms of 
component failure, repair, and system reconfiguration. There are a 
multitude of models which can be used for component failure and repair, 
and reconfiguration. Two good references to the mechanics of reliability 
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theory are [Shooman, 1] and [Green and Bourne, 2]. 


^ Consider a device which begins operation at time 0 and can experi- 
ence catastrophic (i.e. , instantaneous) failure to a non-operational 
state. Let the probability of failure of this device occuring in the 
interval [0,t] be 

F(t) = prob. of failure in [0,t] (1.2.1) 

This is the definition of the failure distribution function [Shooman, 1]. 
Define the hazard rate as 


z{t) = 


dF(t) 

dt 


( 1 . 2 . 2 ) 


1 - P(t) 

from [Shooman, 1]. The hazard rate is the incremental failure probabil- 
ity at time t, given that the device is operational at time t. Now, 
suppose the hazard rate of the device is independent of time; i.e., the 
probability that the device will fail sometime in a time interval ' 
starting at the present time is independent of how long the device has 
been operational. This constant hazard rate 

z(t) = c (1.2.3) 

results in the exponential failure distribution shown in Figure 1.1. 

The constant hazard rate is a close approximation to the actual hazard 
rate of many devices. For example, the transistor has a hazard rate 
similar to that shown in Figure 1.2. This type of function is quite 
common [Shooman, 1]. Early failures in Region I of Figure 1.2 are 
failures during the ”burning-in" of the device; they are associated with 
poor assembly, defective materials and other random fluctuations in the 
manufacturing process. Failures in Region III are due to the wearing out 
of elements in the part. Region II is relatively constant and closely 
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approximates the constant hazard rate function. In a large system, parts 
are generally "burned- in" before assembly is completed; therefore, the 
system begins operation in Region II. As the system ages, periodic 
maintenance removes old parts before the hazard rate rises in Region III. 
Therefore, the assumption of a constant hazard rate is usually justified. 
In this report, the constant hazard rate function is used exclusively. 

This xs due not only to its broad applicability, but also to the fact that 
any non-constant hazard rate requires a reliable control system to keep 
track of the starting times of the system's mode of operation. 

In the discrete-time case , to which this report is confined exclu- 
sively, the hazard rate becomes the probability of failure (or repair or 
reconfiguration) between time t and time t+1. For a system with many 
operating modes , the probability of being in a given mode at a given 
time, given some past probability vector over the various operating 
modes, can be modeled by a Markov chain. If is a vector 

7T^ e (1.2.4) 


where there are L+1 operating inodes, then is propogated in time by 


I^ttl = ^^t 


(1.2.5) 


where 


P = (p,. . ) e R 


L+l X L+1 


( 1 . 2 . 6 ) 


and 



prob. of system being in mode i at time t+1, given it 
was in mode j at time t 

(1.2.7) 


(see [Paz, 33). 
the hazard rate. 


The probability p^^ is the discrete-time equivalent of 
and is time-invariant. In the future, a time-invari- 


ant Markov chain will be assumed as a model of the modes of operation 
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and the statistics of the random switchings between modes . 

It is now necessary to define precisely these modes of operation 
and their dynamic transitions . The terms system conf iguratibn and 
system structure will be used. 

Definition . 4 : System Structure ; A possible mode of operation for a 

given system; the components, their interconnections, and the informa- 
tion flow in the system ^ ^ given time . 

Definition 5 ; System Configuration ; The original design of the system, 
accounting for all modeled modes of operation, and the Markov chain 
governing the configuration , or structural , dynamics (transitions among 
the various structures) . 


An example of three possible structures for a given system is shown 
graphically in Figure 1.3. In this report, structures are referenced by 
convention by the set of non-negative integers 

I = {0,1,2,3, . . . ,l} ‘ (1,2,8) 


The configuration for the design illustrated in Figure 1.3 is depicted 
graphically in Figure 1.4. The nodes of the graph in Figure 1,4 
represent the system structures of Figure 1.3. The edges of the graph 
represent probabilities of transfer from one node to another, and are 
elements of the matrix P. 


p. . = prob. structure i at time t+1 given structure j at 

1+1, D 1 

(1.2.9) 


The state of the system configuration at time t is the structure in 
which the system is operating at that time. 
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Figure 1.1: Exponential failure distribution. 
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Figure 1.2; Typical hazard rate function for a transistor. 
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Figure 1.3: Three hypothetical system structures. 
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k(t) = structural state at time t . (1.2.10) 

k(t) e I (1.2.11) 

This structural state evolves in time to form the structural trajectory 
(of length T+1) 


= (k(0),k(.l), . . . ,k(T)) (1.2.12) 

In general, this structural trajectory is a random variable with apriori 
probability of occurance 

= \(0) ,0*^k(l)k(o')^k(2)k(l) *” ^k(T)k(T-l) (1.2.13) 

(Figure 1.5). 


1.2.2 Control Theory . 

t 

In this report, only linear systems with a quadratic cost index 
are considered. At this time, any more general formulation is of dubious 
value in that the linear quadratic problems can demonstrate many of the 
fundamental concepts of reliable control system design. It is 
doubtful that any other formulation could be solved without the knowledge 
gained from the linear quadratic solutions presented in the remainder of 
this report. As a further restriction, perfect observation of the system 
state is assvuned. The general class of linear systems discussed in 
this report is of the form 


-t+1 ~ -k(t) -t ■^-k(t)-t 

The set of pairs (A, /B.,, ) describe the possible system structures, 

~~ K K 


where 


k(t) e I 


(1.2.15) 


The remainder of the configuration is specified by the Markov chain 
equation (1.2.5). The objective of this research rs to develop control 
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laws which account for the possible structural trajectories (1.2.12) 
while minimizing some function of the cost. The cost function for a 


T-1 


given random state and control trajectory ( ^ ) is 


T T T 

x^Qx_^ +U^Ru^ + X Qx 
— t— — t — t 1. — T— — T 


(1.2.16) 


The function of the cost which is minimized is generally taken to be the 
expected value of over all possible structural trajectories x^. It 
is shown that this class of optimization problems yields solutions 
which are sensitive to both system performance and system reliability, 
as modeled in the configuration. 

In the remainder of the report , ' only variations in the B^matrix , 
or actuators are considered. An actuator is a device which transfers 
the control input to the system dynamics. The actuator in the ^-matrix 
may model a physical linkage, such as is found on the control surfaces of 
aircraft, or, for example, the effectiveness of a tax reduction on the 
economy.. A single actuator may fail in many different modes. For 
example, the ^-matrix can be of the form 

2-0 = [^0 ^ -1 I ‘ ‘ * I -J ^ (1.2.17) 

where the 's are actuators which may fail to an actuator having zero 

gain with a .failure probaility per unit time p^; 

b^ £ (1.2.18) 

Then the system structures representing modes of failure would be modeled 
as &-matrices having at least one zero column. 

This class of linear models can also be used as a model for self- 
reorganizing systems; the only restriction is that the reorganization, 
or reconfiguration, process must be modeled with a constant hazard rate. 
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An important aspect of this research is the study of various types 
of redundancy. At present/ the effect of redundancy on system, performance 
is poorly understood. There are two basic types of redundancy; component 
redundancy and functional redundancy. Component redundancy is the- use 
of two or more identical components (in this report/ actuators) for the 
same task. A good example is provided by equation (1.2.17). Suppose 
two actuators, b . and b . , are identical. If b . fails (Equation. (1.2.18)), 

—1 —j —1 

b j is still operational, and vice-versa. In order to lose the function 

of actuators ^ . and b . , both actuators must fail; this event will 
1 D 

have a lower probability of occurance than the event of the failure 

of b^ ; kj were not in the configuration the function of actuator 

b . would be lost. 

— 1 ■ 

The problem W'ith component redundancy in control theory is how 
should the allocation of control resources be allocated to the redun- 
dant components, and how should the component reliabilities affect the 
choice of an optimal control law? The control methodologies presented 
in this report answer the question for a specific class system confi- 
gurations . 

Functional redundancy implies the overlapping of function of two 
or more components in a system. If one of the components fails, part 
of its function is still performed by the other (redttndant) component ( s ) . 
Functionally redundant actuators are modeled in this report in the same 
way as component redundancy. The functional redundancy is accounted for 
in the expectaion of the cost index over the structural trajectories. 

The dynamics of repair and reconfiguration are all modeled in this 
report as exponential failure distributions (constant hazard rates) . 
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AS an example , if two actuators (b ^ 

and b ^ ) are in a system configura- 

tion and can 
respectively 

each fail with probability p and p per unit time, 

^0 ^1 

, to an actuator with zero gain (0) , then the configuration 

dynamics 

are 

, assuming independence of failures: 


-0 

= 

'^0 

I^ii 


(1.2.19) 

•B>i 

= 

[0 1 

-1 ^ 


(1.2.20) 

^2 

= 


1 03 


(1.2.21) 

^3 


[0 1 

01. 


(1.2.22) 

-0 

-j- 

^1 

with probability 

p (l-p„ ) per unit time 
^0 ^1 

(1.2.23) 

-0 

-> 

^2 

with probability 

p (1-P- ) per unit time 
^1 ^0 

(1.2.24) 

-0 


^3 

with probability 

p p per unit time 

^1 ^2 

(1.2.25) 

^1 

-j- 

-3 

with probability 

p per unit time 

(1.2.26) 


-5“ 


with probability 

p per unit time 

(1.2.27) 


From this information, the Markov chain transition matrix can be formed 

0 0 0 


P = 


1-P^ -Pf +P^. 

0 1 ^0 1 


p (1 - P- ) 1-P 

0 1 2 


0 


P. (1-Pf ) 
1 0 




'0 1 


0 1-P^ 


Pf Pf 
2 1 


0 

0 

1 


(1.2.28) 


Repai-r is considered to be component replacement, and is modeled in the 
same manner; e.g., 

^ B Q with probability p^ . 


(1.2.29) 
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Reconfiguration is the restructuring of the (actuator) configuration to 
compensate for failure, and is modeled as 

— l"^ —4 probability (1,2.30) 

where is a new actuator configuration which will be used on reconfi- 
guration after failure. 

The methodologies presented allow the study of the effects of 
failure, repair and reconfiguration on the optimal control of linear 
systems; they yield a quantitative analysis of the effectiveness of a 
given system design, where effectiveness is a quantity' relating both 
the performance and the reliability of a configuration design (see 
Appendix 1) . 

1.2.3 General Nature of Results. 

There are three classes of reliable controller methodologies: 

I) Passive (Robust) Controller Design 

II) Active (Switching) Controller, Passive Configuration Design 

III) Active Controller, Active Configuration Design 
This report concentrates entirely on classes I) and II) . Class III) 
methodologies are much more difficult to study. The Markov chain models 
of configuration dynamics which work in classes I) and II) do not hold 
in class III) ; as yet, there is no satisfactory way to model the 
configuration dynamics of a system in such a way that the control rules 
are well-defined. 

Class I) methodologies are passive designs. These designs account 
for the occurance of failures in the initial selection of the control 
law; on-line, this class of designs does not use any current 'estimate of 
the structural state of the configuration. The design is "consejrvative" 
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in that it continues to stabilize the system without regard to the current 
structural state , A special sub-class of these designs is the robust 
controller designs. A robust controller will stabilize any structure of 
the system without regard to the configuration dynamics ; i.e.-, if the 
system remains in any structural state forever, it will still be 
stabilized. The class I) methodologies are represented by the 
non- switching gain methodology of Chapter 5. 

Class II) methodologies are active control lers in some sense, 
they are adaptive. Prom knowledge of the system's past, these controllers 
switch their control law on-line in order to compensate for what they 
estimate to be the correct structural state. For deterministic systems, 
these controllers can be determined analytically. For stochastic 
systems, the optimization problems cannot be solved analytically in 
general due to the dual control effect [Fel'dbaum, 4-7]. Thus, 
suboptimal control strategies must be used. The class II) methodologies 
are represented by the switching gain methodology in Chapter 3 and 
its suboptimal extensions in Chapter 4. 
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1.3 Relations with Previous Literature . 

This research is based on a background knowledge in both reliability 
theory and systems theory. Both mathematics and probability theory are 
fundamental in these fields. As general references to the techniques 
used in this report, in real analysis, and measure and integration 
theory, [Rudin,8] , [Segal & Kunze, 9], and [Halmos,10] are good; in 
matrix theory, [Gantmacher ,11] is the standard reference. In probabil- - 
ity theory, [Bauer, 12] and [Doob,13] are definitive; expansions bn the 
theory of Markov chains are found in [Chung, 14] and [Derman,15] . 

There are several good texts on reliability theory; of these, 

[Greene & Bourne, 2] and [Shooman, 1] are possibly the best. [Cox, 16] 
and [Corcoran, 17] demonstrate the current methods of the scheduling and 
use of redundancy in reliability technology. Other good treatments are 
found in [Barlow and Proschan,18] and [ Gnedenko, 19 ] . 

In control theory, a good treatment of. the deterministic linear 
quadratic regulator problem is found in the IEEE Transactions Special 
Issue edited by [Athans,20], and in [Athans s Falb,21] . The dual 
control problem is described in [Fel'dbaum, 4-7] and several other 
publications . 

Previously, several authors have studied the optimal control of 
systems with randomly varying structure. Most notable among these is 
[Wonham,22] , where the solution to the continuous time linear regulator 
problem with randomly jumping parameters is developed. This solution is 
similar to the discrete time switching gain solution presented in 
Chapter 3. The random parameters are restricted to be a continuous 
time Markov chain. The' most notable difference is that in [Wonham,22], 



30 


the assumption is made that the controller has perfect information about 
the present state of the random process on-line . Tlie solution switches 
gains in a linear state feedback control law whenever the (Markovian) 
random parameter jumps. In the discrete time switching gain solution 
presented in Chapter 3, the control law is determined from pash observa- 
tions -which allow the deduction of the exact state of the random para- 
meter process, and then the random parameter may switch values according 
to the statistics given by the Markov chain. Thus, the control may be 
applied to one of a number of possible structures at the next time 
instant. In Wonham' s development, the optimal control law is matched 
specifically to one structure. The analogous continuous time version 
to the- switching gain solution of Chapter 3 would be to assume on-line 
perfect observation of the random parameter with a fixed time delay. 
Wonham 's result has no such time delay. 

Wonham also proves an existence result for the steady-state optimal 
solution to the control of systems with randomly varying structure. 

This result is based on conditions of stabilizability of each system 
structure and observability of each structure with respect to the 
cost functional. The conclusion is only sufficient; it is not necessary 
for existence of a steady-state solution. Similar results were obtained 
in [Beard, 23] for the existence of a stabilizing gain, where the 
structures were of a highly specific form; these results were necessary 
and sufficient algebraic conditions, but cannot be readily generalized 
to less specific classes of problems. 

The time-varying solution of [Wonham, 22] is computed using a set of 
coupled Riccati-like matrix equations. The coupling is in the form of 
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a linear term in the solution to the matrix equations added to the normal 
linear quadratic Riccati equation. The solution can be precomputed by 
solving the coupled Riccati-like equations off-line; the control law is 
then switched on-line to a gain which corresponds to the current state 
of the Markov process. The optimal solution requires perfect knowledge 
of the structure. 

In reality, the structure is seldom known perfectly, and a noisy 
observation of the random process leads to a dual control problem. 
Although much of Chapter 3 is based on the fact that the controller can 
obtain the structural state with one-step delay in the deterministic 
discrete time problem, this report makes the connection, for the first 
time, of the existence of a steady-state switching gain controller with 
that system's reliability and effectiveness. 

[Sworder, 24] has developed, using- a version of the stochastic 
maximimi principle, an optimal feedback control law for a class of' linear 
systems with jump parameters which is almost identical to that of 
|Wcxiham,22]; the coupled Riccati-like equations are identical except for 
notation. The only difference is Sworder 's assumption that the random 
process is instantaneously observable from a set of sensors which are 
unaffected by the choice of the control law. Using this assumption, 
Sworder avoids the problems of dual control. 

Sworder also comments on the usefulness of linear system models 
with jiimp parameters in modeling possible failures in the system 
[Sworder ,24 ] . [Ratner & Luenberger ,25 ] derive a control law for a 
continuous time linear system. The system has one failure mode , and a 
maximum number of renewals (repairs) can take place. The objective is 
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to determine apriori the optimal time intervals in, which the system 
should operate in the failure mode, and the optimal - control law, given 
the mode of operation, over a finite time interval. The failure process 
is assumed to have an exponential failure distribution (constant hazard 
rate) ;■ the renewal process is controlled, and is not random. The 
control law is of the switching gain type, and the solution is in the 
form of two coupled Riccati-like matrix equations quite similar to those 
in [Wonham,22] and [Sworder,24] . The optimal control policy and the 
optimal renewal policy can both be calculated off-line. This class of 
problems is further investigated by [Sworder,26] to determine over what 
region immediate renewal is the optimal policy. Both of these papers 
illustrate examples of class III) control methodologies ; the structural 
state as well as the system state is under the influence of the control- 
ler. The simple structure of the class of systems studied by [Ratner & 
Luenberger ,25] allows a solution. There is need for much more work in 
this area. 

Still a third approach to the problems associated with multiple- 
structure systems is given in [Bar-Shalom & Sivan,27]. Here, the 
measurements of the system state are corrupted by additive noise. The 
open-loop controller and the open-loop feedback controller are derived 
using dynamic programming. Knowledge of the present state of the random 
process governing the system configuration is not assumed. Therefore, 
the (optimal) closed-loop controller would be a dual control law. The 
open- loop controller assumes no on-line measurements of the system state; 
the open-loop feedback controller assumes future on-line measurements 
and thereby improves its performance. There is little correlation 
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between this paper and the research on which this report is based. 

[Winner, 28] developed a suboptimal control scheme, which allowed 
for imperfect observation of the random parameter process, known as 
multiple-model, adaptive control. In this method, the parameters, could 
only take a discrete set of values, a cause of recent disfavor, as MMAC 
does not always work well when the parameters vary continuously and are 
approximated by the mathematics. Similar work has been done in [Pierce s 
Sworder,29] . The MMAC methodology is optimal one step backward from the 
final time, as is the switching gain methodology in the example of 
Chapter 2 when applied to systems with additive white control noise. 

The dual problem of state estimation with a system with random 
parameter variations over a finite set was studied in [Chang S Athans,30]. 
It is shown there that the optimal estimator consists of a geometrically 
increasing set of Kalman filters, one for each possible structural 
trajectory of length t+1 at time t, and an averaging process to compute 
the minimiim mean-square error estimate from the filter estimates. It 

is also shown that when the parameter process is Markovian, a bank of 

2 . . . 

N estimators is optimal, where there are N possible values of the 
parameters. Each estimator is then conditioned on the possible values 
of the parameters at the two previous time instants. 

Recently, the robustness of the linear quadratic regulator has been 
studied in depth. This work is described in [Wong, et . al. ,31] and 
in [Safonov & Athens, 32] . A long-standing problem with the linear 
quadratic design methodology has been the lack of analogs to the various 
stability and robustness criteria of classical systems theory. This 
research was aimed at characterizations of robust solutions to. 
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specif icall'y, the linear quadratic regulator. Supporting research is 
reported in [Safonov & Athens, 33], [Wong & Athans,34] , [Wong, 35], and 
[Safonov, 36] . The research in this report is related to the robust 
controller problem, but the approach is different in that the performance 
criteridn- is modified to account for possible variations in structure, 
such as those caused by failures, rather than depending on certain 
properties of the linear quadratic regulator solution to guarantee 
robustness. In this research, the concept of stability is related to 
the existence of a finite cost solution to the non-switching gain 
problem. For a specific class of configurations, this approach solves 
the robust controller problem (Chapter 5, Section 9) . 

The existence of an uncertainty threshold for the non-switching 
controller of Chapter 5 , that limit on parameter uncertainty beyond 
which no controller can stabilize the system, is proven for an one- 
dimensional example. This work is similar to the work by [Athens, 
et . al . ,37 } on the Uncertainty Threshold Principle and the related 
papers by [Ku & Athans,38] and [Ku, at. al . ,39 ] . This research is 
reported in Chapter 2, Section 7. 

Lastly, parts of this research have been presented in an unpub- 
lished form at the 1977 Joint Automatic Control Conference in San 
Francisco, and piiblished for the 1977 IEEE Conference on Decision and 
Control Theory in New Orleans [Birdwell & Athens ,40]. 
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1 . 4 Summary of Main Contributions . 

There are two major contributions of this research. First, the 
classification of a system design as reliable or unreliable , for the 
deterministic variable actuator linear system in Chapter 3, has been 
equated with the existence of a steady- state switching gain and cost 
for that design. If the steady-state switching gain does not exist, 
then the system design cannot be stabilized ; hence, it is unreliable . 

The only recourse in such a case is to use more reliable components 
and/or more redundancy- Reliability of a system design can therefore 
be determined by a test for convergence of the set of coupled Riccati- 
like equations (3.3.6) as the final time goes to infinity, 

A similar result holds for the non— switching gain methodology of 
Chapter 5. Here, the system design is classified as reliable^ or 
unreliable with respect to a constant gain linear feedback control law, 
depending on the convergence, or divergence, respectively, of equation 
(5.6.16) as the final time goes to infinity. If equation (5.6.16) 
converges to a limit cycle, then that limit cycle produces a stabilizing 
cyclic steady-state gain. 

The second major contribution lies in the robustness implications 
of the non-switching gain methodology. Precisely, a constant gain for 
a linear feedback control law for a set of linear systems is said to 
be robust if that gain stabilizes each linear system individually , i . e . , 
without regard to the configuration dynamics. The problem of determining 
when such a gain exists, and of finding a robust gain, can be formulated 
in the context of the non-switching gain methodology. As a result, the 
non-switching gain methodology gives an algorithm for determining a 
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robust gain for a set of linear systems which is optimal with respect to 
a quadratic cost criterion. If the algorithm does not converge, then 
no robust gain exists. 

The following Section of this Chapter will outline the remainder 
of this report. 

1.5 Outline of Report . 

In Chapter 2, several one— dimensional examples are examined as 
a clarification and motivation for the methodologies presented in 
Chapters 3 through 5. In addition. Chapter 2, Section 7, deals with 
the relationship between the Uncertainty Threshold Principle and the 
existence of a steady-state solution to the non-switching gain problem. 

Chapter 3 develops the optimal solution to the class of problems 
described in Section 2 of this Chapter. The solution is labeled the 
switching gain solution because the gain of a linear feedback control 
law switches in response to the exact observation of the system 
structure with one-step delay. 

Since Chapter 3 deals entirely with deterministic systems , and the 
switching gain solution does not extend optimally to the stochastic 
case , Chapter 4 presents some suboptimal methods which can be used to 
extend the switching gain solution to stochastic problems. Two 
methodologies are presented- One (hypothesis testing) is based entirely 
on estimation of the structure. The second (dual identification) uses 
the dual effect of the control law to determine more precisely what the 
structure is with the next observation. The optimal control law would 
have some characteristics of both methodologies, as is shown by example 
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in Chapter 2, Section 5. 

Chapter 5 derives a control law which ignores any on-line info 2 nna“ 
tion which might be gathered about the structural st-ate, and results 
in a non-switching gain solution used in a linear feedback control -law. 
The stability of this non-switching solution is explored, along with 
the existence of a steady-state solution, in Secion 7. In Section 9, 
the robustness issue is addressed, and the non-switching methodology is 
used to define an algorithm which can determine the existence of a 
robust gain and calculate an optimal robust gain with respect to a 
quadratic cost functional, when one exists. 

Chapter 6 focuses on the issues of computer-aided design and the 
application of the non-switching gain methodology to design problems. 
Two examples are used to demonstrate the effectiveness of the non- 
switching methodology in design. 

Chapter 7 reviews the results described in the report and suggests 


new directions for future research. 
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CHAPTER 2 

CLARIFICATION AND MOTIVATION OF RESEARCH 
2.1 Introduction . 

The purpose of this Chapter is to motivate all subsequent more 
general Chapters with simple one-dimensional examples. In particular, 
in Section 2, a one-dimensional problem is formulated and solved to 
illustrate the optimal (switching gain) deterministic control for 
linear quadratic systems with variable actuator configurations. 

The effects of process noise on this solution are examined in 
Section 3. The dual effects which occur in the stochastic systems 
motivate the suboptimal approaches described in Chapter 4. 

The possibility of steady-state control of variable actuator 
configuration systems with a single linear independent control law 
is discussed in Section 6, motivating the work on the non-switching 
gain solution and robust control laws in Chapter 5. In addition, 
the possibility of existence of a steady-state stabilizing linear 
feedback control law with constant gain is compared with the work on 
the Uncertainty Threshold Principle [Athans,et.al. , 37] in Section 7. 
Section 7 contains the only case of this report where exact algebraic 
conditions for the existence of a steady-state solution have been 
derived. Unfortunately, these results do not readily extend in an 
analytical manner to higher dimensions. 

The question of existence of a steady-state solution to these 
problems is of great importance. A system design is defined to be 
reliable with respect to a certain class of control laws if there 
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exists a control law from that class for which the infinite time 
cost incurred using that control law is finite. Since the switching 
and non-switching, gain solutions are- the- optimal solutions for- their 
respective classes of control laws, if they incur an infinite cost, so 
will any other control law from that class. In addition, since the 
switching gain solution is the optimal control law for the determin- 
istic problem, a system design is termed deterministically reliable , 
reliable if and only if the incurred infinite time expected cost 
is finite. 

In the next Section, a one-dimensional example is presented 
which will' be used to motivate the remainder of this report by 
examining the ramifications of the switching and non-switching gain 
solutions through their specific application to the example. 

2.2 A Simple Example — The Optimal Solution . 

The following one-dimensional example is used to demonstrate the 
switching gain methodology presented in Chapter 3, and to show that 
the general stochastic problem is analytically intractable. All proofs 
and derivations are given in Appendix 2. 


2.2.1 Problem Statement. 

Let the discrete-time system be one-dimensional with one control 


variable u^ and state variable x^ related by 

X = ax, + b, u 
t+1 t k t 


( 2 . 2 . 1 ) 


The value of the control multiplier (b^) is a random variable which 


takes on one of two discrete values at each time t. 
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( b if k = 0 
b^(t) = j 


( 2 . 2 . 2 ) 


1/b if k = 1 

The random process k(t) is governed by the Markov chain represented 
by 

IT = PTT (2.2.3) 

— 1+1 1 

where 

^ • (2.2.4) 


IT ^ e R 
— t 


P = 


^^11 

^12 

.^21 

^ 22 . 


(2.2.5) 


At any given time t, the following sequence of events occurs: 

I) is observed exactly, computed, and k(t-l)is 

set to 0 or 1 depending on where k(t-l)is the 

variable representing the Markov chain; 


II) b, may change values to b, ; 

k(t-i) k(t) 


III) is applied. 


J = ^ (qx^ + ru^) 


( 2 . 2 . 6 ) 


For any given sample path, the performance index is given by 

TP 

s 

where {o,1,...,t} is the time set over which the system is to be 
controlled. The objective of the control problem is to minimize the 
expected cost-to-go at time t, given by 

V(x^,k(t-1) ,u^,t) = (qx? + ru“) |k(t-l)l (2.2.7) 


where the expectation is taken over all possible sample paths of 


k(T) , t < T < T. 
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2.2.2 Summary of Solution. 

From Appendix 2.1, we find that the optimal control is given 
by 


u.. 


•"o.t 

+ "o,t ”l,t 


( 2 . 2 . 8 ) 


where 


2^t 


TT 


TT 


0, t 

1, t 


= P TT 

--t-1 


(2.2.9) 


Thus, the control law is linear in the state x^, and switches between 

two precomputable gains, depending on the value of k(t-l). 

Given X. , , and , 

t t-1 t-1 


^t-1 


^t"^t-l 
if — — = b 


u 


t-1 


^t-“^t-l 


( 2 . 2 , 10 ) 


u 


= 1/b 


t-1 


and k(t~l) = 0 if ^ — t-l~^*^ U ’ - 

The optimal cost-to-go is 

V*(x ,k=i,t) = x^S. (2.2.11) 

where S and S are propagated backward in time by the following 

O # t _L / "t 


equations: 
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Assuming k=0 at time t, then TT ^ 


=o,t = "3 + 


[r+p^j^b^So^t+i+Pji iW) Sj, t^.j.1 ^ 


+ Pll . a ^ ^o,t+l 




■^^211®“ " 2 ,^2, I "l,t+l (2.2,12) 

b[r+p^^b SQ^^^jL"^p 23 _S^^^^^/b ] 


Assuming k=l at time t, then tt ^ = [p ^2 ^22^ ' 

,2 


=i,t - * 


‘Pl2°^^0.t-l-l+P22 
tr+Pi2b^So ,.^1+P22 (1/b^) Si_ ^ 




^ , Pl2*‘’^O.CH.l'^P22<°'^’h.t.|.l . , 

^22 1^ ^ v.2^ ‘ C /k2i I 1/t+l (2.2.13) 

Note from equation (2.2.8) that u^ switches from one linear gain 
to another, depending on the value of — thus,, this solution depends 
on an exact knowledge of x^. If knowledge of x^ is corrupted by measure- 
ment noise (or, if u^ is corrupted by control noise) , then it will be 


shown by example that this becomes a dual control problem. 
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2.3 The Dual Control Effect. 


+ i*, / ^ u , + Co. 
t+1 t k(t) t t 


To demonstrate the difficulties encountered when white process 
noise is present, -the optimal solution for the one dimensional 
example is derived over the time interval {o,l,2} with additive white 
control noise present. The system is now represented by 

(2.3.1) 

C. is discrete time white noise with zero mean, ] = H6 , 

t ^t^T t-T 

probability distribution p(^), and is uncorrelated with and k(r) 
for T_< t. 

* * 

Thus, the problem is to find u^ and such that the expected 
cost-to-go is minimized. 

From Appendix 2.2, the optimal control one step back in time 
(at t=l) is 


Ui = - 


2 T^idlDb^ 


, qa 


r t 


Li= 


7T^(l|l)bJ 


(2.3.2} 


where ir^(ljl) is the probability that k^ = i, given the information 
set expected, this control is of the same 

form as is the deterministic control law, equation (2-2.8), since 
there is no benefit in trying to determine k^ more accurately through 
the use of a special control value. In other words, there is no dual 
control effect at t = T^-1 (in this example, t=l) . 

At t-0, the situation is different. Now, the optimal control will 
force the system to supply more information through the state at t=l 
than it normally would in the absence of the process white noise 



44 


In order to compute u^, a numerical minimization of a numerical 

* 

integration (in general) must be performed. Thus, u^ is the 


solution of 


V (Xq, 0) = min 


! ^ . 

uixii ( X g + u 


,q + u^r + :iq 


ko=o 


X? 1' 

1 */R(x^) 


q(l+a^) 






r + 


i 



V 

dp(xJki,k„,Zo)p^^^^ 

O 

O 


Li=0 


(2.3.3) 


where 


ir^dll) 



p(Xi-aXp-b.Ug)Tr.^O 

i P(X;j_-axQ-b.uQ)7r.^Q 




(2.3.4) 


and p (Xj^ jk^,k^,ZQ) is the probability measure of x^ over R(x^)., 
the range of x^, given k^, k^, and Z^. 

Equation (2.3.3) is very difficult to solve numerically, and 
for any realistically-sized problem would be economically infeasible. 
For the limited amount of computation that has been done with equation 
(2.3.3), the dual control effect is evident from Table 2.1. Note 
that as the process noise variance increases, the trend is for the 

•k 

control Uq to increase. This is due to the need for a larger control 
to lessen the effect of noise on future estimations of the structure. 
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Table 2.1 


The optimal control versus and H. 



Uo -( H- 3 ) 

(5=6) 

o 

( — t 

11 

:iii 

c O 
b 

-2.0 

2.3170089 

2.3188635 

2.3201611 

-1.6 


1.8550055 

1.8559061 

-1.2 

1.3898305 

1.3907551 

1.3912676 

-0.8 

0.9255912 

0.9259997 

0,9261950 

-0.4 

0.4606236 

0.4606920 

0.4607206 

0.0 

-0.005 

-0.005 

-0.005 

0.4 

-0.4706236 

-0.4706920 

-0.4707206 

0.8 

-0.9355912 

-0.9359997 

-0.9361950 

1.2 

-1.3998305 

-1,4007551 

-1.4012676 

1.6 

-1.8635511 

-1.8650055 

-1.8659061 


* - calculation did not converge due to numerical errors 


The system used in the calculations is described by equation 
(2.2.1) where 
a = 2. 

k(t) is 0 or 1 

bo = 2. 

bi = .5 
q = 3. 
r = 1. 



Table 2.1 is only . intended to demonstrate the difference in the 
optimal control laws at time 0 for a two-stage process; numerical 
accuracy is not assured. Specifically, the values of -.005 for 
Uq (Xq = 0) are highly doiibtful, as well as the consistent 


asymmetry between positive and negative values in the Table. 
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2.3.1 A Special Case. 

It is interesting that for one specialized probability distribu- 
tion p(^), when the optimal control is large enough, the optimal 
solution is identical with the deterministic solution of Section 2. 

Prom Appendix 2.3, assuming 

^ 2 ^ ' for 

P(0 = i (2.3.5) 

I 0 otherwise 

* 

as shown in Figure 2.T, if from the deterministic solution (equation 
2.2.8) satisfies 

I (bj^ ~ > 2 /^ for k’^ i (2.3.6) 

* ^ 

then Uq is also the solution to the stochastic control problem. 

Physically, because the noise is amplitude limited, it is easy 
to exactly deduce the structure if the control is large enough. 




L: A probability di, 

white noise. 
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2 . 4 Existence of a Steady-State Solution . 

Although, as will be stated in Chapter 3, little can be said 
about the existence of a steady-state solution to the general n-dim- 
ensional switching gain problem, for the one-dimensional example, 
exact conditions for the existence of a steady-state solution can 
be found. They are in the form of two simultaneous algebraic equations 
which can be solved analytically. 


r - Pii I « - 


b[p^^ab+P2j^{a/b)hI \ 2 

Pl^b^+P2i(l/l»^)b 


+ P 


21 


a - 


p^^ab+P2j^(a/b)h 

>’[Pll‘=^+P2lh/b^> 


(2.4.1) 


hF = p 


12 






+ p 


22 


a - 


p^ 2 ^b+P 22 (a/b)h 
b[pi2b^+P22bA>^3 


(2.4.2) 


The equations are derived in Appendix 2.4. In these equations the 
variables F and h are defined as 

r - Urn (2.4.3) 

t->--<» 0 , t+1 

and 

^1 t 

h = lim — - (2.4.4) 

°0,t 

whenever both S and S increase without bound as t as defined 

in equations (2.2.12) and 2.2.13). 
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S 


0,t 


Since F is the limiting value of the ratio of the next value of 
to the present value S , it is necessary that 

U f t T J, 


r > 1 


,(2.4.5) 


for 


S a> 

0,t 


(2.4.6) 


Similari'ly, if S has a limit, then F can have a maximum value of 
Oft 

1. Therefore, a test can be made on the solution (h,D to equations’ 
(2.4.1) and (2.4.2) for the existence of a steady-state solution; 

If 


h 7 ^ 0 or 0° (2.4.7) 


then 


%,t' ^ ^ ^ 

(2.4.8) 

S , . converge if F < 1 

Oft J- / 1 

(2.4.9) 


and there is no conclusion if F = 1. 

By way of eliminating all possibilities, as an aside, a limit cycle 
to the solution of equations (2.2.12) and (2.2.13) cannot occur by 
Lemma 1 of Chapter 3. 


2.5 Conclusions on the Switching Gain Methodology . 

The purpose of the last three Sections on the one-dimensional 
switching gain example was to clarify the approach of this phase of the 
research, and to motivate the approach of Chapters 3 and 4. In this 
Section, some implications of the one -dimensional example will be 


discussed. 
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2.5.1 Implications of the Dual Control Effect. 

It was shown in Section 2 that the optimal solution to the deter- 
ministic class of variable actuator linear quadratic control problems', 
i.e., the switching gain solution, is conceptually straightforward, 
although computationally complex off-line. Unfortunately, in Section 3, 
it was demonstrated that the optimal solution of the stochastic version 
of the same problem is infeasible. (Witness the problems of calculating 
the two-step optimal solution.) Therefore, since the switching gain 
deterministic solution is essentially the only solution which can be 
described analytically, the research involved in developing the 
n-dimensional switching gain solution is justified. This is exactly 
what is presented in Chapter 3. 

It then remains to investigate any extensions (which will of 
necessity be suboptimal) which may be made to the switching gain 
solution to adapt the solution to the stochastic problem. In Chapter 
4, a start is made in that direction. These are two basic routes 
to follow: The various hypothesis testing algorithms in combination 

with the switching gain solution, and a formulation developed in 
Chapter 4 which gives the control vector a dual effect; the control 
is changed to increase the accuracy of the estimation algorithm. 

The optimal control would use techniques from both categories, as the 
dual effect is clearly seen in Table 2.1. 

2.5.2 Existence of a Steady-State. 

Although for the one-dimensional example, it is possible to 
determine the condition for convergence of the Riccati-like equations 
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(2.2.12) and (2.2.13), this method does not extend to the n-diiuen- 
sional solution. It is at present unknown under what conditions 
the Riccati-like equations for the n-dimensional problem converge; 
therefore, there is little comment on conditions for convergence in • 
the remainder of this report. 
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2,6 A Simple Example — The Non-Switching Solution . 

In the prevxous sections of this Chapter, motivation was given 
for the development of the optimal (switching) solution to the linear 
quadratic ■ variable actuator configuration control problem. 

Several problems with the method were pointed out in Section 5. 
Specifically, the methodology does not extend optimally to the stochas- 
tic case due to the dual control effect. Secondly, the increase in 
on-line complexity over the usual linear quadratic control problem 
is significant, especially in the suboptimal stochastic schemes. 

In many instances, a stabilizing solution to this class of 
control problems is desired which exhibits the same complexity as 
does the usual linear quadratic controller. For instance, it may be 
desired that a control law stabilize a system without requiring 
error detection strategies and switching to a new form upon detection 
of failure. A subclass of these problems occur when a robust gain 
(one which stabilizes each configuration without regard for the 
dynamics of structural changes) for a set of linear systems is 
desired. The first problem within this subclass deals with the 
existence of such a gain. The second problem deals with the choice 
of an optimum robust gain with respect to some cost index. 

In the following Subsections, an example of non-switching gain 
methodology is given as an illustration of the concepts; since the 
derivations are quite complex, proofs are deferred until Chapter 5, 
where the entire development of the non-switching solution is presented 
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The following formulation is only for the steady-state solution; 
in Section 7, the conditions for existence of the steady-state solu- 
tion will be given and related to the Uncertainty Tihreshold Principle 
[Athans,et. al. , 3 


2.6.1 Problem Statement 

In Chapter 5, the non-switching control problem is solved for 
linear systems with variable actuator configurations and quadratic 
cost. It was stated in the conclusion of the previous Section that 
a relationship exists between the existence of a steady-state solution 
and the Uncertainty Threshold Principle, In this Subsection, the 
existence of a steady-state non-switching solution to the one dimen- 
sional example presented in Section 2 will be studied to illustrate 
this relationship. 

The system to be used is 


X ^ , = a x_^ + b, u^ 
t+1 t k t 


( 2 . 6 . 1 ) 


where x, a, b^ and u are scalars, k can be either 0 or 1, and t takes 
on integer values. 

b if k=0 




( 2 . 6 . 2 ) 


1/b if k=l 

The index k represents the structural state of the system, and 


is a random variable with statistics generated by the Markov chain 
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(2.6.3) 


P = 


P 1-P 

1-p p 


(2.6.4) 


where ]T . . is the probability that the structural state is i at time t, 

2. fX- ' 

given some initial condition ir (T , . ) . 

— init 

The inf inite- time , or steady-state non-switching control problem 


is formulated by specif ing that the solution u ^ is to minimize the 


cost of a trajectory f ^t ^ t=T given by the sum 

°° init 

qxt + rut 

^ ^ init 


(2.6.5) 


2.6.2 Summary of Solution 

The solution is computed, from Chapter 5, equations (5.7.17) 


and (5.7.18), when it exists, as the solution (S ,S,) of 

o 1 


- 2 / /. 

= - ;rr2— TT- 


+ 


(bSQ+S^/b)^(r+b^S^) 


5 (b^'S^+S^/b"') +r 4 (^ (b^S^+S^/b^) +r) ^ 


+ (1-p) ( - 


+ q 


(bS|j+s^A3) 


(bSQ+s^/b)^(r+S^/b^) 


(%(b^S^+S^/b^)+r)b 4(Js(b^S^+S^/b^)+r)^ 


( 2 . 6 . 6 ) 


2 1 ,. (bSQi-S,/b)bS^ 

Si = M 1 Sq- ~TT— T- 


(bS^+S^/b)^(r+b^SQ) 


ig (b“SQ+S^/b"-) +r 4 ih (b^S^+S^/b^) +r) ^ 


+ P(S^- 2 


'bSo+S/b)S_^ 


(bS^+S^/b)^ (r+S^/b^) 


(%(b‘’S^+S^/b^)+r)b 4(%(b^S^+S^/b^)+r) ^ 


+ q 


(2.6.7) 



and the control is given by 
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* 



(bSg+Si/b)a 

(r+%(b^S^+S^/b^)) 


( 2 . 6 . 8 ) 


Note that the steady-state solution is a linear feedback control 
law with a constant gain which is pre-computable using equations 
(2.6.6) and (2.6.7). The on-line implementation of this solution has 
the same complexity as does the usual linear quadratic steady- state 
solution. 
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2 . 7 Existence of ^ Steady-State Solution and the Uncertainty 

Threshold Principle . 

In this Section, the existence of a steady-state solution to 
equations (2.6.6) and (2.6.7) is related ho the Uncertainty Threshold 
Principle [Athans et . al . ,37 ] . This Principle states that for a 
certain class of systems, there exists a threshold, or bound, on the 
degree of uncertainty in the system dynamics beyond which no control 
law will stabilize the system. Furthermore, it is noted in 
[Athans et . al . ,37 ] that there does exist a "minimizing" control even 
though the infinite- time cost in infinite. 

For the non-switching gain class of controllers, it will be 
shown in this Section that, at least for the one-dimensional example 
of Sections 2 and 6,- such a threshold does exist; furthermore, it will 
be explicitly calculated. In addition, it will be demonstrated that 
the non-switching control gain -converges even when no finite cost 
steady-state solution exists. 

2.7.1 Formulation of Existence Problem. 

The question is now asked: When does the steady state solution 

exist? I.e., when is the cost, given by 

J = 3g(SQ + S^)Xq (2.7.1) 

finite? 

This problem is solved by showing when the solution does not 
exist. 

Allowing 

Sq ^ ~ (2.7.2) 


and setting 
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h = lim 


l,t 


tr^-oo 0 , t 

S> 


r = lim 


0,t 


(2.7.3) 


(2.7.4) 




where S and S are the values of the r.h.s. of equations (2.6.6) 
U / J- / 1 

and (2.6.7) iterated backwards t times from an initial value S. =0/ 

1,0 ^ 


equations (2.6.6) and (2.6.7) become 


r = a ^ 




(b+h/b)b 

%(b^+h/b^) 


(b+h/b) 

,, 2 , ,^ 2 , 2 
(b +h/b ) 


+ (1-p) h 1 - 


(b+h/b) 


hVo^+h/b^)h 


(b+h/b) ^/b^ 

2 2 2 
(b +h/b ) 


(2.7.5) 


hr - a I (1-p) 1 - 


(b+h/b) b 

^ (b^+h/b^) 


(b+h/b) ^b^ 

2 2 2 
(b +h/b ) 


+ P 


H'-i 


(b+h/b) 


®i<b^+h/b^) b 


(b+h/b) ^/b^ 

2 2 2 
(b+h/b ) 


(2.7.6) 


2.7.2 Summary of Solution. 

Equations (2.7.5) and (2.7.6) have 5 solutions. The solutions of h 
and r of interest are: 


For p / 

h = -(p(b^(6-2W)-3b®-3) + ((2b^-2)p-b'^+l)V 

+ (4b®-2b'^+2)p^+b®-2b‘^+l)/( (2b'^+2)p^-2pW) (2.7.7) 

r = a^(-p[b^(2p^+4p-2)+(b®+l) (p^-2p+l)]^+(b^+l)p^) 


/((b^+l)^(2p-l) 


(2.7.8) 
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where 

V = [b^(p(4-4[b‘^(2p^+4p-2)+b^(p^-2p+l)+p^-2p+l]'^) 

+2p^-2) +b® (5p^-2p+l) +p^-2p+l] (2.7.9) 

and 

W= [{b®+2b^+l}p^+{~2b®+4b^-2)p+b®-2b'^+l]^ (2.7.10) 

For p = h: 

h = 1 (2.7.11) 


r = 


2 2 
(b^-1) 

2(b'*+l) 


2 


(2.7.12) 
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2.7.3 Graphical Illustration of Solution. 

Equations (2.7.7) through (2.7.12) are tod complex for much 

information to be gleaned from study. Therefore, their significance 

is demonstrated graphically in this section. 

These equations are used to compute the absolute values of a 

versus b and p above which no stabilizing non-switching control exists; 

i.e., since F is the limiting ratio of . to S. , ,, what threshold 

O , t 0 , t+X 

value of I a} yields F = 1? Since the system (2.6.1) is a discrete 
time one, this threshold quantifies how unstable the open-loop system 
must be for there to be no stabilizing solution. This quantity is 
called the uncertainty threshold value of { a | . For the case p = %, 


^threshold ®^sy to compute from equation (2-7.12) 


' threshold 


[2(b^+l)]^ 


For p ^ 



(b^+1) r (2p-l) 

/ (p ( (b^+l)p- [b^ (2p^+4p-2)+ (b®+l) <p^-2p+l]^) ) 


(2.7.14) 


A plot of a 


threshold 


versus p (long axis) and b is shown in Figure 2.2. 
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Figure 2.2: a 


threshold 


versus p, b. 
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The ln(b) axis is used because af , , ^ , is symmetric with respect 

to ln(b) around zero ( I a | (b) = I )- b varies 

from e to e p varies from p = 1 to p = -JDI . Note that 

la| , , , ^ “ as b 1 and/or p 0. This is because as b ->■ 1, the 
‘ ' threshold 

system looks more and more like 

x_^ ^ = ax^ + bu^ (2.7.15) 

t+1 t t 

which is controllable for all values of a. As p ^ 0, the system is* 

switching more and more rapidly between the two structures; therefore, 

each structure has less time to influence the system unfavorably and 

the system becomes easier to control, leading to a , , , 

' ' threshold 


2.7.4 Best Control with Infinite Cost. 

Although the cost may be- infinite, a finite gain control exists. 
Prom equation (2.6.8), and allowing and ^ h, the control 
becomes 


,^2 ^ „ 2, "t 
(b +h/b ) 


(2.7.16) 


Note that the control gain does not depend on q or r, but only on p, 
a and b, as in the work with the Uncertainty Threshold Principle. A 
plot of h versus p (long axis) and b is given in Figures 2.3a and 2.3b, 
in the same manner as for F. Note that as p ■*> 0^, h “ (except at 
b = 1) . For this boundary, we rely on a symmetric argument, switching 
the roles of and S^, since we only know that ->* 

An interesting symmetry exists in h with respect to p. If h is 
defined as 


h 


lim h 
b->0 


(2.7.17) 



76265AW025 


62 



DBIGINAL 
0E POOR 


Figure 2.3a: h versus p, b. 







Figure 2,3b: h versus p, b. 
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then 

h = (2.7.18) 

P 

Letting p = % + x, 

<2.7.19) 

and 

h(x) = ^ — (2.7.20)' 

h(-x) 

Thus, ln[h(p)] is symmetric around p = .5.. This solves the boundary 

problem, because as p-^ 1, h ^ 0 (except at b = 1) , and the condition 

Sq -»• <» is satisfied (S^^ ->-0) . Since h is symmetric, and h(p,b) -J-h(p) 

for p->"0, the solution is well-defined at p = 0. 

In Figure 2.4, the control gain divided by a, g,‘ is plotted as a 

function of p and b. 

* 

u^ = -gax^ _ (2.7.21) 

'f' *p- 

Note that as p 0 (and h ^ g -»■ b, and as p ->■ 1 (and h 0) , 
g -»■ 1/b, and that bQ = b and b^^ = 1/b. Thus, as p -s- 0^, the optimal gain 
tends towards the deadbeat controller for the system in structural 
state 1, and as p 1 , the optimal gain tends towards the deadbeat 
controller for the system in structural state 0. 
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Figure 2.4; g versus p, b. 


ORIGINAL PAGE IS 
OF POOR QUALOT 



66 


2.7.5 Conclusion. 

In this Section, the steady-state properties of the non-switching 
solution to a specific example of actuator failure were studied, and 


were related to the Uncertainty Threshold Principle. In particular, 
the existence of an uncertainty threshold has been established, and 
with the help of the high degree of symmetry in the example, the values 
for I ^ 1 ' given b and p, were calculated. It was also shown 

that the best control with infinite cost is a function only of a, b and 
p, a situation analogous to the solution obtained in the papers on the 


Uncertainty Threshold Principle [Athans et. al.,37]. 


An analogous solution to that presented here should exist for the 
switching gain problem, and in fact, the rudiments of such a solution 
are given in Section 4. As a guide for future research, it would be 
interesting to compare the two methodologies on the basis of these 
solutions. Unfortunately, it is mathematically intractable to extend 
this result to the multivariable case, although another approach may 
be found. 


2.8 Summary . 

The unifying issue in this research is the interrelationship 
between the issues of control and reliability. Section 7 brushes on 
the question of when a system design is considered a reliable design. 
In Chapter 3, a reliable design will be defined as one in which the 
steady-state switching gain solution exists. Therefore, questions 
concerning the existence of such solutions become quite important. 
Unfortunately, little headway has been made in the development of any 
simple test for the existence of the steady-state solution. Only in 
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Section 7, in the specific case of the non-switching gain solution, 
for a specific (relatively trivial) example, and in Section 4 for the 
same example with the optimal solution., have conditions for existence 
of a steady-state been resolved. In section 7, these conditions are 
given explicitly; in Section 4, they are given as the solution to two 
simultaneous equations. For the general n-dimensional problems in the 
remainder of this report, existence can only be tested by iteration of 
the solution equations. 
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CHAPTER 3 

THE SWITCHING GAIN SOLUTION 

3.1 Introduction . 

In this Chapter, a control methodology for linear systems with 
quadratic cost criteria and variable actuator configurations will be 
developed which accounts' for the failure, repair and reconfiguration 
of the actuators by switching the control gain on detection of a 
change in configuration. This problem is viewed as a control problem 
rather than as the traditional estimation problem. Therefore, a 
deterministic model is assumed, except for the random changes in 
configuration, which are modeled by a Markov chain. This methodology 
has the advantage that all gain and expected cost calculations are 
^one off-line . The gains switch on-line with changes in the configura- 
tion, which are observable with one-step delay for almost all values 
of u^ (i.e., except for a set of measure zero). In addition, the 
method is useful in the stochastic case, though not optimal, in 
conjunction with identification methods such as hypothesis testing 
and dual identification, which will be described in Chapter 4. The 
gain and expected cost calculations can be used as an evaluation 
technique in computer-aided design of linear systems. An example 
would be in trade-off studies of various redundancy configurations 
with respect to performance, reliability, and system effectiveness. 

The disadvantages of the technique as it is presented here are that it 
requires perfect measurement of the state and that only multiple 
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actuator configurations are considered. The multiple sensor configura- 
tion problem should be dual to this work. Changes in the A matrix 
are a minor extension; however, the general, problem allowing variations 
in both the actuators and the observers would be a major result. 

Previously, several authors have studied the optimal control of 
systems with randomly varying structure. Most notable among these is 
[Wonham,22], where he develops a solution to the linear regulator 
problem with randomly jumping parameters in continuous time. The 
solution assumes apriori that the controller has perfect information 
about the present state of. the random parameter process. Little work 
was done on the steady-state existence problem. 

The solution presented in this Chapter is analogous to that of 
Wonham's; however, the discrete time formulation of the problem allows 
the controller to observe exactly with one step delay the value of the 
Markov parameter process. Thus, it is shown that for the discrete- 
time process, the optimal controller is not dual. 

In addition to this conclusion, this research makes the connection, 
for the first time, of control and system reliability and effectiveness. 
This is the unifying concept in the entire report, and has been discuss- 
ed in detail in Chapter 1. 

The procedure for determining the existence of a steady-state 
solution to the switching gain control problem divides system designs 
into two classes; If a design allows a steady-state solution, then 
that solution is stabilizing (see Section 7, Chapter 5); therefore, 
that design is classified as a reliable design . On the other hand, if 
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no steady— state solution exists, then that design is classified as 
inherently unreliable. 

Although no easy test exists for the existence of a steady-state 
solution, the computer can always be used to iterate equation (3.3.6) 
backward in time and check for stability. Therefore, this methodology 
yields a classification of systems into those which are inherently 
reliable and those which are not. 


3 . 2 Mathematical Formulation . 

In this Section, the n-dimensional extension to the one-dimension- 
al switching gain result presented in Chapter 2 will be developed. 

The only, non-trivial task is to prove that the system structure is 


observable for almost all values of the control. The system model is 


“t+1 


(3.2.1) 

where 



it ^ 


(3.2.2) 

iit ^ 

R™ 

(3.2.3) 

A e 

^nxn 

(3.2.4) 


and, for each k, an element of an indexing set I 
k £ I = {0,1,2, . . . ,L> 


where 


(3.2.5) 

(3.2.6) 


B, £ {b ^ (3.2.7) 

— k —1 r e I 

The index k(t) is a random variable taking values in I which is 


governed by a Markov chain and 
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•nr ^ , = p IT ^ 
— t+l — t 


— t ^ 


L+l 


(3.2.8) 

(3.2.9) 


where tt. ^ is the probability of k(t) = i, given no on-line,information 
r ,.t 

about k(t),, and IT^ is the initial distribution over I. 

It is assumed that the following sequence of events occurs at 
each time t: 

1) is observed exactly 

2) then B , , . switches to 

k(t-l) — k(t) 

3) then is applied. 


The control interval is assumed to be 
{0/l,2, . . . ,T> 

and the cost function is selected as 


T-1 

T T T 

r=u 


T 


(3.2.10) 


(3.2.11) 


The objective is to choose a, feedback control law, which may 
depend on any past information about x ^ or , mapping into u^ 

^ 71 in 

: R — >-R (3.2.12) 

•k 

<t> . i — ►U. (3.2.13) 

— t — t — t 

such that the expected value of the cost function from equation 
(3.2.11) 

= E [ I TT^ ] (3.2.14) 

is minimized over all possible mappings ^^at . 
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3.3 The Switching Gain Solution . 

Normally, a control law of the form (3.2.13) must provide both 
a control and an estimation function in this type of problem; hence 
the label dual control is used. Here, the structure of the problem 
allows the exact determination of k (t-1) from ^ and ^ 

for almost all values of u ^ ^ . This result is stated and proved in 
the following theorem. 


Theorem 1: For the set {b . } . where the B. 's are distinct, the 

— k kei — k 

set {x has distinct members for almost all 

values of u ^ . 

Proof: See Appendix 3.1. 


Ignoring the set of controls of measure zero for which the 
members of 

^-k,t+l^k=0 (3.3.1) 


are not distinct, then for (almost) any control which the optimal 

algorithm selects, the resulting state compared with the 

members of the set (3.3.1) for an exact match (of which there is only 

one with probability 1) , and k(t) is identified as the generator of 

that matching member x , . . 

— k,t+l 

Since perfect identification is the best any algorithm can achieve, 

Jr Jr 

the optimal control law u^ = ^^(2£^) can be calculated with the 
assumption that k(t-l) is known, since this is the case with probability 
one. Thus, this solution will be labeled the switching gain solution , 
since, for each time t, L+1 optimal solutions are calculated apriori, 


and one solution is chosen on-line for each time t, based on the past 
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measurements x ^ — t ^ which yield perfect knowledge of 

k(t-l) . 

Dynamic programming will be used to derive the optimal switching 
gain solution. At each time t, the expected cost-to-go using the 


control sequence 


-t ' -t+1 ' -t+2 ' • 


/ u 


T-1 


(3. 3,-2) 


and given the value of k(t-l) is defined as 


V(x ^ ,u^ ,k(t-l) ,t) 
T T 


+ (x ,k(t) ,t+l) |k(t-l)} 


(3.3.3) 

where * denotes the optimum value and v’' is the optimal value of V. 
Then, by dynamic programming 


V (3£ ^ fk (t-1) , t) = min 




' 7? ^ T ^ 


+ '''(M-t+l) I (3.3.4) 


It is proved, from Appendix 3.2, that 
v*(x^,k(t-l),t) = 


(3.3.5) 


where the ^ are determined by a set of L+1 coupled Riccati-like 
equations (one for each possible configuration) : 


:,t = j ^ p s . 


t+i 


^%k-i,t+l-i — ^ ^ik -i— i,t+l — i 

T 1 I 

2^ik^i^i,ttl ^-*-2 


-1 


(3.3.6) 
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The optimal control, given k(t-l) - k, is 


- k,t 


R + 


* P-i B.S. 

A-f( ik —1— i,t+l 1 

1=0 


^ p bT S . . . , B . 
4-i ^ik — 1— i,t+l— 1 
1=0 


-1 


Writing 

* 


- k,t 




then 


-k,t 


= - j-R + ^ p b'^. S . B . 
— ^ik — 1 ~l,t+l — 1 


(3.3.7) 


(3.3.8) 



^ik -i-i,t+l- 


(3.3.9) 


Thus, is a switching gain linear control law which 

depends on k(t-l). The variable k(t-l) is determined by 

k(t-l) = i iff x^= Ax ^ ,+ B . , (3.3.10) 

— t t“l — 1— t~l 

Note that the ^ , ' s and the optimal gains can be computed 

X f iC f *C 

off-line and stored. Then, at each time t, the proper gain is selected 
on-line from k(t-l), using equation (3.3.10), as in Figure 3.1. 


3.4 Discussion of Results . 

The solution in section 3 is quite complex relative to the struc- 
ture of the usual linear quadratic solution. Each of the Riccati-like 
equations (3.3.6) involves the same complexity as the Riccati equation 
for the linear quadratic solution. In addition, there is the on-line 
complexity arising from the implementation of gain scheduling. In 
Chapter 5, a non-switching gain solution will be presented which has 
an identical on-line structure to that of the linear quadratic 
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solution, but has similar of£-line computational complexity to that of 
the switching gain solution. Depending on the system requirements, 
either solution could be used; the non-switching gain solution is 
suboptimal, but requires less on-line complexity. This trade-off may 
favor the non-switching solution in some cases. 

A steady-state solution to equation (3.3.6) may exist, but the 
conditions for its existence are unknown. The steady-state solution 
would have the advantage that a time-invariant set of gains result. 
Thus, only one set of gains need be stored on-line, instead' of requir- 
ing a set of gains to be stored for each time t. Since the steady- 
state solution is simply the value to which equation (3.3.6) converges 
as it is iterated backward in time, at present, the equations can 
be iterated numerically until either they converge or meet some test 
of non-convergence. Unlike the non-switching solution presented in 
Chapter 5, the possibility of limit cycle solutions in the switching 
gain computations is excluded by the following lemma; 

Lemma 1 : If the optimal expected cost-to-go at time t is bounded 

for all t, then equation (3.3.6) converges. 

Proof ; See Appendix 3.3. 

Once again, it is stressed that the existence of a steady-state 
solution to the switching gain problem establishes a division of 
system designs into those which are inherently reliable and those 
which are unreliable . Even though conditiois to test for the exis- 
tence of the steady-state solution are unavailable , software can be 


used with iteration for the test. 



In Section 5 , some numerical examples are given to illustrate 


the switching gain solution. 
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3.5 Examples. 

In this Section, a two-dimensional example is presented with three 
different switching gain solutions to illustrate ‘the switching gain 
computational methodology. The computer routines which are used in 
the calculation of the switching gain solution are listed in the 
Appendix. The primary subroutine is READY; it calls WEIGHT. Any other 
routines which are used are from the standard ESL siibroutine library. 

The main program RDYMAIN is used to call READY. 

Example 3.1 is a two-dimensional system with four structural 
states corresponding to the failure modes of two actuators. In this 
example, failure of an actuator is modeled as an actuator gain of 
zero. Thus, the four structures are: I) Both actuators working (B ^ ) ; 

II) One actuator failed (^^and ), and III) Both actuators failed 
(^3 ). The system is controllable in all structures except for the 
sturcture represented by ^3 . 

Actuator failures and repairs are assumed to be independent events 
with probabilities of failure and repair, per unit time, of p^ and p^, 
respectively, for both actuators. 

In Example 3.1, the matrixes Q and R are the quadratic weighting 
matrices for the state and the control u^ , respectively. The 
matrix ^ is the Markov transition matrix, which is calculated from knowl- 
edge of the system configuration dynamics, represented graphically 


in Figure 3.2. 
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Figure 3.2: Markov transition probabilities for Example 3.1 
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There are three Cases to Example 3.1. Each Case assumes a different 
failure rate and repair rate for the actuators. Case i) has a high 
probability of failure and a low probability of repair, relative to 
Cases ii) and iii) . The switching gain solution is not convergent for 
Case i) ; the gains themselves converge, but the expected costs do not. 
Only configuration state 0 is stabilized with its corresponding gain. 



Cases ii) and iii) both assume more reliable actuators than does 
Case i) . Both Cases ii) and iii) have convergent switching gain 
solutions. Therefore, both Cases ii)' and iii) represent reliable 
configuration designs , while Case i) is unreliable . This difference 
is due entirely to the different component reliabilities. Equivalently, 
Cases ii) and iii) are stabilized by the switching gain solution, while 
Case i) is not. Note that in this Example, s tab iliz ability is not 
equivalent to stability in each configuration state, or robustness ♦ 

For this example,, no robust gain exists because the system is 
uncontrollable from configuration state 3. 

Cases ix) and iii) are also presented in Chapter 5, where their 
non- switching gain solutions are given. According to the theory, it 
should be more difficult to stabilize a given system with the non-switch- 
ing gain than it is with the switching gain, because of the optimality 
of the switching gain solution. This is demonstrated for this example; 
in Chapter 5, the non-switching gain solution to Case ii) is not 


convergent. 



Example 3.1: 


A = 


2.71828 0.0 

0.0 .36788 


' 1.71828 

1.71828' 

II 

«T 

o 

f 

o 

■ 

1.71828' 

-.63212 

.63212_ 


o 

• 

o 

.63212. 


'1.71828 

r 

O 

O 

-3 

o 

0 

1 

o.o' 

.-.63212 

O.OJ 


[0,0 

0.0. 



'14. 

8.' 


'l.O 

0.0 ' 

= 



R = 




8, 

6. 


O 

• 

o 

1.0 , 


l-2p^+P£ 

(i-Pf)Pr 

(i-Pj)Pr 

2 

Pr 


Pf Cl-Pf ) ^ 

i-Pf-p^+PjPr 



-Pr) 

Pf{l-Pf)" 

PjjPf 

l-Pf-P,-^PfP^ 

Pr<l- 

-P^) 

2 

Pf 

(l-Pj,)Pf 

(l-P^IP, 

l-2p^+P^ 


The system dynamics are 


it+l 


2-t 




k(t) £ {0,1,2,3> 

The cost, which is to be minimized, is 


J = E 


GO 

T I 

2£t 22it I 1 

t=0 
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Example 3.1, Case i) 

_ 


— — 

Pf = .3 

.49 


^0 

, 

II 

« 

.21 


^1 

IT = 


— 



.21 


^2 


.09 




Non-Convergent ; but 

gains converge: 


1 

-.9636 

1.094 X lO"^ 


II 

o 

-.9134 

-5.835 X 10“^^ 



-.9234 

1.740 X 10~^ 



-.8699 

-5.136 X lo"^ 



-.8094 

.9186 X 10“^ 


^ 2 = . 

- 1.020 

-4.0,5 X lO"® _ 



-.9636 

.7353 X lO”^ 



-.9134 

-3.923 X 10“^ 

- 

Stability: 


Configuration 

Stable 



o 

ml 

o 

yes 



1 (B ^ ) 

no 



2 (B^) 

no 



3 (B 3 ) 

no 
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Example 3.1, Case li) 


Pf = -1, 


P.. = -9 


.81 


^0 

.09 


^1 

.09 



.01 


^3 





Convergent Coupled Riccati Equations: 


G . = 
— 1 


-.8390 

-.7752 


.04222 

-.09914 


S . = 
— 1 


25.57 

8-611 


8.611 

6.398 


for i = 0,1, 2, 3 


Stability : 


Configuration Stable 


0 (B Q ) ■ yes 

1 (B ^ ) no 

2 (B 2 ) no 

3 (B 3 ) 


* no 
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Example 3.1, Case iii) 

p, = .01, p = . 98 
f r 


.9799 


— 

■^0 

.009999 


■^1 

.009999 


■^2 

.0001020 


■^3 



_ 


Convergent Coupled Riccati Equations : 
-.7558 .1270 ■ 

-.8073 -.1786 



15.88 

8.1051 

^0 = 

8.105 

6.137J 


-.7060 

.1186 

^1 = 

-.8441 

-1.723 


16.06 

8.074 


8.074 

6.143_ 


-.8375 

.1090 


-.7543 

-.1669_ 


'16.31 

8 . 199 ' 

^2 = 

8.199 

6.158 



85 


^ 3 = 


-.7863 
7926 


.1023 

-.1619 


16.54 

8.170 

4 


Stability: 


8.170 
6. 162 


Configuration Stable 

0 Q ) yes 

1 ) no 

2 (B 2 ) no 

3 {^2 ) no 
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3*6 Summary , 

In this chapter, the optimal solution to the linear control 
problem with variable actuator configuration was developed. It was 
shown that the optimal solution uses a linear switching feedback gain 
which depends on the previous configuration. This configuration is 
directly computable from the past measurements; this fact allows the 
development of the switching gain solution by eliminating dual con- 
trol considerations. The exact measurement of the configuration with 
one-step delay holds only for the deterministic case, where there is 
no corruption of the state or control observations by noise. 

In Chapter 4, the use of the switching gain methods- will be 
demonstrated for stochastic problems in conjunction with two different 
forms of identification! Hypothesis testing and dual identification, 
a technique for "pushing" the control variable out of the noisy 
region, when the noise is amplitude limited, to obtain an exact 
identification of the system structure. 
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CHAPTER 4 

EXTENSIONS TO THE STOCHASTIC CASE 
4.1 Introduction . 

In Chapter 3, the optimal solution to the deterministic linear 
quadratic control problem with variable actuator configuration was 
developed. It was also demonstrated that the optimal solution of 
the general stochastic linear quadratic problem- is hopelessly complex 
in Chapter 2. Therefore, in this Chapter, extensions to the deter- 
ministic solution to allow its operation in a stochastic environment 
will be studied. 

From the derivation of the switching gain solution, whenever 
the structure of the system is known perfectly with one step delay, 
and if it is assumed that it will be measured perfectly at the next 
time instant, the optimal solution is the deterministic switching 
gain solution. In designing a suboptimal control system, a method 
of identifying the system structure is used, with the assumption that 
the identification is perfect, and the appropriate deterministic 
gain is selected. 

Two conceptually different methods of structure identification 
will be presented in this Chapter. The first is classical hypothesis 
testing. It is the easiest to implement, although extensions to 
n-step hypothesis testing can be made which are very complex. The 
second method is labeled dual identification ; the expression is used 
because it takes advantage of the dual effect of the control law to 
guarantee perfect identification. In this method, a perturbation 
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(which may or may not be that small) to the deterministic control is 
introduced which separates the effect of amplitude limited white 
control noise from that of the system structure. As a worst case 
control law, this perturbation would be applied at each time instant, 
but in practice, it would only be applied once every n time instances 
so that its overall effect on system performance would be lessened- 
In the next Section, the system model will be described, and the 
hypothesis testing identification algorithm will be presented. 


4 . 2 Hypothesis Testing Identification ♦ 

The system model used here is the same as in Chapter 3, but with 
the exception that additive white noise is introduced into the 
dynamics ; 

For the hypothesis testing identification method, is assumed to be 

zero mean white noise with probability distribution p (Q . It is 

assimed to be uncorrelated with k (t) and x . Perfect measurement of 

— t 

the state is retained. 


The basic hypothesis testing method is very simple: At each time 

t, one of L+1 hypotheses is chosen, where each hypothesis is 

: k(t-l) = i (4-2.2) 

With each hypothesis H^, there is a probability of being 
correct, given the measurement and the past information: IT (t—1 1 1—1) , 

the probability distribution of k(t-l), given the measurements through 
, Then the updated' probability (see Appendix 2) Tr^(t-ljt), the 
probability of k(t-l) = i, given all measurements through , is 
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given by 




P(2Lt 



- )lT^{t-l|t-l) 

“ -j ^t-1 


(4.2.3) 


Hypothesis is assumed to be correct if 
7T^(t-l|t) > TT^(t-llt) for all j i 


(4.2.4) 


Ties are resolved arbitrarily. Then, given the correct hypothesis H^, 
the corresponding deterministic optimal switching gain is used to 
compute the control at time t 

u=G.x (4.2.5) 

“t “i,t-t 

as in equations (3.3.8) and (3.3.9). 

The probability distribution is then propagated with the Markov 
chain equation 


TT(t|t) =^Ti^(t-l|t) (4.2.6) 

and the process repeats. 

This algorithm can work well if there are significant differences 
in the effect of the control variable between configurations. ’When 
the differences are slight, a mistracking will result until the errors 
are large enough to be detected through equation (4.2,3), The method 
does not exploit any of the dual effect of the control variable on 
the measurement of the configuration. The method presented next does 
use the dual effect to identify the correct structure. Analytically, 
it cannot be said which method is best, as the optimal control law 
will lie somewhere between the two. It is possible to extend the 
hypothesis testing procedure to n-step hypothesis testing where a 
hypothesis is made about the last n values of k(t) and is then tested. 
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since this investigation is not within the primary scope of this 
research/ it is left as an open problem for future research. It is 
also possible -that a- combination of hypothesis testing and dual identi- 
fication may be used to gain some of the advantages of both methods; 
dual identification yields fast identification of the correct structure/ 
while hypothesis testing does not sacrifice control of the system 
while there is a high probability that the structure is correctly 
identified. 


4.3 Dual Identification . 

The underlying concept of dual identification is to periodically 
change the control in order to increase the accuracy of identification 


of the structure- In the limiting case, the control is changed 
enough to guarantee perfect identification of the current structure 
with the next observation. For this case only amplitude limited noise 
is considered. The system model is 

it+i = +Mi^ (4.3.1) 

where is ^-dimensional white noise which takes on values in the 
unit sphere with distribution p(_§^) and is uncorrelated with and 

k(t) . M is an nx£ matrix which defines the ellipsoid in which 
contains 

Normally, if no identification were to be performed, and if k(t-l) 


were known, the optimal deteanninistic switching gain G . . from 

equation (3*3.9) would be used to compute u* . 


* 



— k(t-l) ,t-t 


(4.3.2) 


In dual identification, the goal is to compute a gain offset u , ^ 

1 / t 



such that when the control 


Et = (4.3.3) 

is applied to the system, identification of the structure k(t) with 

the observation x is guaranteed. To accomplish this, note that, 
t+1 

for a given B, , , will be in a bounded convex set determined by 

— k —t+1 

B and M. Thus, 

1C 


-t+1 


Ax =B u + Mg (4.3.4) 

t JC ^ ~ t 

jj, 

and can be any element in the unit sphere S (R ) . Therefore, 

perfect identification of k(t) is guaranteed if no two of the domains 

of corresponding of the B^'s have a non-empty open intersection 

That is, the following condition must be satisfied for each pair of 

I 

B , ' s and every 5 and 5 o S (R ) ; 

Jc , j- z 

X 


'it ^ 0 


(4.3.5) 


This condition is the same as 

II (Bj^^ ~ )U^ II >2 


if (1.V - )u ^ ^ N.(M) 

kf k^ t 


otherwise , 

(B - B, )u, 0 (4.3.6) 

^1 ^2 ^ 

where M is the generalized inverse of M and N (M) is the nullspace 

of M. Note that the inequality of (4.3.6) can be relaxed to equality, 

since the intersection of the two domains of x , , would only be at 

— t+1 

the point of tangency, a set of measure zero in either domain. 

The objective is to choose such that (4.3.6) is satisfied 

for all pairs B and B , in the reachable subset of all actuator 



configurations. The reachable subset refers to the subset of configu- 


rations which have a non-zero probability of occurance at time t, 
given that the configuration was t-1. This is the same as 

the condition that 


B . is in the reachable subset from B , , ^ , 

— 1 — k(t-l) 


(4.3.7) 


Suppose that there are J configurations in the reachable subset from 

B , . . Then there are J<J+l)/2 pairs of configurations for which 

condition (4.3.6) must be satisfied. Also, since u affects the 

JL f 

state X , it is reasonable to minimize its effect. Therefore, 

— 1+1 

since the effect of u is modified by B . . , it is reasonable to 

1 / u K {X.) 

minimize the norm of u . Thus, the minimization problem is formu- 

1ft 

lated siabject to the constraints (4.3.6) • 

2 

min 


in|| I 

-l,t 

subject to 


4 - <0 


where 

D. = M (B. - B.) 

-i+(3-l)*J - -1 -2 

Formulating this as a nonlinear programming problem, the 
Hamiltonian is 


,X) - II u 


^ \(4- ||Dj^ -l,t' 


A > 0 

k — 


Aj^ = 0 if 4- 

Differentiating H with respect to A_, and solving for ^ 

dt 

function of u^ and the parameter X, 


1/t 


(4.3.8) 


(4.3.9) 


(4.3.10) 

(4.3.11) 


as a 
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M. 

9u 


= 0 = 


1/t 


-i-t - ? 




or. 


u 


l,t 


= II - E V 


T -1 

^ k^k-k ' 
k k 


E T * 

A D D u 
k-k-k-t 


(4.3.13) 


Now, using (4.3.13) in the constraint equation (4,3.11) 

1- ii£^ti+ [I- II" 


(4.3.14) 


Noting that ■ 

U - XU + Xd'^D] “^0*^21] = U + Ad'^D^]"^ (4.3.15) 

then (4.3.14) simplifies to 

-1 - II E, II - E \e!e^ rh* II < 0 ( 4 . 3 . 16 ) 

k 

and if (4.3.16) is a strict inequality, then X^^ = 0. In general, 
a niimerical algorithm must be used to solve for ^ in the set of 
equations (4.3.16); this can be a major drawback to the application 
of this methodology if the on-line computer resources are unavailable. 
Although the computational burden of this technique is a disadvantage, 
dual identification would most likely be implemented in combination 
with a hypothesis testing algorithm. Dual identification would then 
form a test to be performed on the system after some interval of time 
to ensure that the hypothesis testing algorithm correctly tracked the 
configuration. 
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4.4 Examples . 

In this Section, the one-dimensional example of Chapter 2 , Section 2 
is implemented with additiwe white, noise applied to the control input. 
Three suboptimal control algorithms derived from this Chapter are imple- 
mented: Hypothesis testing, dual identification, and hypothesis 

testing in combination with dual identification, which is utilized every 
fifth time instant. The purpose of this example is to illustrate the 
degrading effect of the dual identification algorithm on the system 
state. 

The principle subroutine used to generate the computer simulations 
of Example 4.1 is SWITCH; it is listed in the Appendix. SWITCH calls 
PIG and UCALC, also in the Appendix; any other routines which are used 
are in the ESL subroutine library. 

The system in Example 4.1 has two structures, represented by the 
matrices (b = 2.) and (1/b = .5); the Markov transition probabili- 

ties are given by the matrix The switching gain solution was calcu- 
lated using the software described in Chapter 3, Section 5. Case i) 
of the Example corresponds to the hypothesis testing methodology described 
in Section 2. The additive white noise was amplitude-limited with zero 
mean and variance H = 1. Case ii) of the example demonstrates the perfor- 
mance degradation due to the exclusive use of dual identification. Note 
that the variation among the values of the state and control are larger 
than in Case i) . The advantage of dual identification is that, for 
amplitude-limited white noise, perfect identification of the system 
structure with one-step delay is guaranteed. In Case iii) , hypothesis 
testing is used four-fifths of the time to partially avoid the degradation 



95 


due to dual identification. The control is more effective in Case iii) 
than in Case ii) ; however, for this example, it is not clear that the 
use of dual identification one-fifth of the time is warranted, since a 
performance degradation of Case iii) over Case i) is still evident in 
this particular simulation. More simulation would have to be carried 
out before the proper ratio of the use of hypothesis testing to the 
use of dual identification could be deteirmined. 



Example 4.1 : 


A = 1.414 


= 2.000 


= .5000 


P = 


Q = 3-000 R = 1.000 

'.7 -3 

.3 .7 

Switching Gain Deterministic Solution: 


Gq = -.7569 


= -1.008 

The system dynamics are 


^'t+i • ®k(t)"t 


k(t) E {0,1} 

The cost function which was minimized is 


J = 


E 


r oo 

E 



where 



Structural transitions are of’ the form 


.3 



When dual identification was employed, the control was set to 

* 

u^ = 1. 25 (sign (u^) ) 

This control was the minimum value required to establish perfect 
identif ication . 
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i) Hypothesis Testing 
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iDDual Identificatim 
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i ii) Hypothesis iesting ond Dual Identification 
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4-5 Siiimnary. 


"In this Chapter, two methods have been proposed to extend the 
-deterministic optimal switching gain solution of Chapter 3 to the 
stochastic case. The two methods represent the two fundamental 
concepts of identification: Estimation and dual control. The 

optimal stochastic control law, if it could be computed, would rely 
on both concepts-, using estimation when the control variable is 
large (and the state is far from the origin) and dual control to 
enhance estimation when the control and state variables are small. 

In the dual identification technic[ue presented here, control is 
sacrificed to obtain an exact observation of 'the structure. Thus, 
the system response would be roughly periodic, with the state being 
driven away from the origin in order to obtain, an accurate estimate 
of the configuration, and decaying back toward zero between identifi- 
cations. In the period when the control is not modified, hypothesis 
testing would be used to track the configuration. 
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CHAPTER 5 

THE NON-SWITCHING GAIN SOLUTION 

5.1 Introduction 

■In the previous two chapters, the switching gain solution was 
developed and studied. In this chapter, attention will be focused 
on obtaining a constant, robust, or non-switching gain which solves 
a: variable actuator configuration linear quadratic control problem, 
with minimum cost for this class- of solutions. ‘It must be stressed 
that thiSi is a suboptimal solution; for the deterministic case. 

Chapter 3 gives the optimal solution. The interest in this chapter 
lies in determining a sequence of ‘gains, for a linear control law, 
which do not switch in response to the detection of a change in system 
•structure. For instance, it may be desirable to- ensure the stability 
of a control system -under certain types of failure without creating 
the- complexity necessary "to- detect those failures and compensate for 
them, as is done in the switching gain solution. 

This class of solutions is related to the overall robustness 
problem where fault- tolerant control systems are desired. Although 
not formulated in this manner, the* research described in this Chapter, 
as in chapter 3, is readily extendable to system with variable system 
matrices as well; i.e., where the- system can be represented as a set 
of possible structures (A^^ B^) over some suitable index, even though 
this class of problems is not as directly related to the underlying 
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reliability theme of this report, 

Non-switching gain solutions to the variable actuator configura- 
tion class of problems oan be obtained i-n- different mathematical ways. 
Problem A of Section 3 is reformulated as a deterministic control 
problem (Problem AE ) , and is solved using the necessary conditions of 
the Matrix Minimum Principle [Athens, 41] in Section 5. Unfortunately 
this approach, although yielding the necessary conditions for an opti- 
mum, does not allow an- analytic solution. Therefore, in Section 6, 
a second problem (Problem B) is formulated and solved using dynamic 
programming . 

Section 7 is by far the most detailed and one of the most impor- 
tant sections of the report, along with Sections 8 and 9. In Section 
7, the concepts of stability and cost- stability are defined and' are 
used to prove an equivalence between the infinite- time versions of 
Problems AE and B. In Subsection 7.6, the steady-state solutions for 
both problems are defined. Unfortunately, nothing in the mathematics 
appears to rule out the possibility of limit cycles in the infinite- 
time solution'; this is discussed in Subsection 7.7. When the constant 
steady-state solutions to the two problems exist, it is proved in 
Section 8 that they are identical. This is a very important result, as 
it allows- the steaey-state solution of a complex two-point boundary 
value problem which is much more tractable. 

In Section 9, it is demonstrated that the general robustness problem 
for linear systems (where one wishes to determine a single stabilizing 
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gain for a set of linear systems) is solved in this framework for the 
class of systems with variable actuator configurations. Examples of 
both the non-switching solution to Problem B and the robustness 
result are given in Section 10, and a chapter . summary in Section 11- 

5.2 Problem Statement . 

The objective of the research described in this Chapter is to 
form a methodology which will be used to compute apriori a gain G 
(either time-varyxng or steady-state) which minimizes the expectation 
of the quadratic performance index over a set of linear systems with 
actuator variation and known transition probabilities of structural 
change (Problem A) . The necessary conditions for minimization are 
given which this optimal gain must satisfy; it is shown that these 
conditions result in a complex two-point boundary value problem. 

A second optimization problem is formulated which is based on 
the restriction to non-learning control laws which are precomputed; 
i.e., it is assumed that the control law cannot benefit from knowledge 
of its past. Although this formulation appears to be much weaker 
than that of Problem A, it is shown in Theorem 2 that if steady-state 
solutions to the two problems exist, then the steady-state solution 
to Problem A is stabilizing (in the sense that the mean square value 
of the trajectory is exponentially bounded) if and only if the steady- 
state solution to Problem B yields a system which is exponentially 
stable- This result is very significant, in that a Corollary to this 
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Theorem solves the problem of finding a robust gain for a set of linear 
systems and yields an explicit procedure for its calculation* 

The last- Theorem (Theorem 3) of the Chapter proves that the steady- 
state solutions to the two optimization problems are identical. This 
implies that not only does the procedure mentioned above determine a 
robust gain if and only if such a gain exists, but also that the steady- 
state gain is optimal with respect to the specified quadratic cost 
criterion. 


5.3 Problem A. 


Consider the system 


X ^ , = Ax B , * u , 

~t+l 1 — k(t)— t 


where 


X ^ e R 
— ■}- 


_m 


u ^ e R 
— t 


(5.3.1) 

(5.3.2) 

(5.3.3) 


k(t) el={o,l,2,'*‘,L} (5.3.4) 

I is an indexing set for the possible actuator structures • 

where 

B, e R (5.3.5) 

~k 

k(t) is a random variable with sufficient statistics given by the 


Markov transition probabilities p^^ , where the matrix 


£ - (Pi3> 


(5.3.6) 

is a stochastic matrix, and the initial probability distribution is 
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(5.3.7) 


Since k(t) is assximed to be a Markov chain, the probability vector 
is propagated in time by 


(5.3.8) 


TT = P 
— t+1 1 

where there is no real-time observation with which to update . 

Consider the structure space {b. } indexed by I. Define the 

JC JCoJ- 

structural trajectory x^to be a sequence of element k(t) in I which 
select a specific structure at time t, 

x^= (k(0), k(l),..., k(T-D) 


(5.3.9) 


The structural trajectory x ^ is a random variable with probability of 
occurance generated from the Markov equation (5.3.8). 


T-1 

= n ^ 


t=0 


k(t) ,t 


(5.3.10) 


where the control interval is 

{0,1,2,. .. ,T-1,T> (5.3.11) 

•for the finite time problem with terminal time T. Then for a given 

T-1 — 

state and control trajectory (5.3.1) and 

T-1 

from a sequence of controls (u^ ^ t-0 ’ cost index is to be the 

standard quadratic cost criterion 




T-1 



(5.3.12) 


The admissible controls are restricted to be of the linear feedback form 




(5.3.13) 


* 


i/e. 



(1 0 ...0) or (0 1 0...0) or ... (0 0...0 1) . 
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where the gain matrix restricted to be a function only of time 

and the initial conditions ; i. e . , cannot depend on • The objective 
is to minimize over the set -of admissible controls the expectation of 
(5.3.12), where the expectation is taken over the set of possible 
structural trajectories 


^ n I 

T 


and the set of initial conditions x^ . 


(5.3.14) 


* * 

Thus, the optimal control law u. = ’ G^x^ should minimize the 

‘ — t — t — t 


cost 


= E[J„iTT ^ ] 
T T' — 0 


= E 


T T T I 

x^Qx^+u,Ru, + x^Qx„Mnr„ 
— 1“— t — t 1 — T— — T'— 0 


(5.3.15) 


over the set of admissible controls. 

Since the structure of u, = G^x^ is fixed, the problem is equiva- 

— t — t — t 

lent to minimizing, in an open-loop sense, the cost function 

T^T, T_T__ , 1^ 

x^Qx^ + x.G.RG.x. +x^Qx_ ir^ 

— t— — t ~t— t 1— t- — T ' — 0 


E[J,pllo3 = E 


(5.3.16) 


with respect to the gain matrix , t=0,l, . . . ,T-1. Equation (5.3.16) 
is simply obtained by substituting equation (5.3.13) into equation 
(5.3.14) . 


5 . 4 The Method of Solution . 

The matrix minimum principle [Athans,41] will be used to determine 

•k 

the necessary conditions for the existence of u^ (or equivalently, 

* 

G ^ ) . To solve the problem using the matrix minimum principle, the 
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formulation presented in the last section must be converted into an 
equivalent deterministic problem. For this purpose, let the- initial 
state Xq be a zero mean random variable which is independent of any 
structure . Let 


Lq = Etx^xJjmQ] = ECXqXqI 


(5.4.1) 


be the convariance matrix of x ^ - 

Defining the covariance of as 

A T I 

Lt = E[x^x^|tLq] 

then, by direct calculation, we obtain 


(5,4.2) 


L L L 

Z = * • • 7 P. . P. . •** P. . TT. „ 

^ i 4^0 i 4^0 x4'o ^t-i^t-2 ^t-2S-3 "-o'° 


t-1 t-2 


t -1 1 r ^ 

n (A+B £ . ) n (A+B G ) 

j=0 j ^ J [j=0 -"j . 


Similarly, if we define 


= E[x 21,- lk(t-l)=i,Ti2^ ] 

X f w w ^ t- ”” W 


then, we deduce that 


(5.4.3) 


(5.4.4) 


t-1 t-2 t-3 


■ t-2 

•TT. (A+B . G ) n (A+B G ) Z 

0 L3=0 3 J 


• n (A+B^ G ) (A+B^G^_^> 
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Remark 1 ; At this stage, an equivalent deterministic problem (Problem AE) 

will be defined with state (2 . . )^ _ for t > 0 and state E „ at t = 0. 

— X , t 1=0 — 0 

The system dynamics are then defined by equations (5.4.6) arid (5.4.7). 

Definition (Problem AE) : For the system with matrix state (2 . , n 

for t>0 and for t = 0 with dynamical equations (5.4.6) and (5.4.7) 
and matrix control , minimize the equivalent deterministic cost 
over 



(5.4.9) 
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Note that since the expectation in equation (5.3.13) is over all 
structural trajectories x and the initial also, 

- j (5.4.10) 

T T 

The symbol will be used exclusively in the future. The one-stage, 
or instantaneous, cost at time t is 

= tr[E^ (Q + G^RG^)] (5.4.11). 

Problem AE is completely deterministic in the state (E . , )^ „ t 

— X , t x=0 — 0 

and control . 

At this point, the minimization will be decomposed into two parts 
using the Principle of Optimality [Athans and Falb,21]. The first 
minimization is over the interval {l,2., . . . ,T-l} , and for this the matrix 
minimum principle will be used. The resulting solution will depend 
in general on the choice of and on the initial conditions and 

* 

Let V (.Gq) he the optimal cost resulting from the use of G^ and 

* * * - ^ 
the optimal sequence » • • • ^ — t- 1 inteirval 1.1,2, — ,Tj- 

The second minimization is then over G ^ of the cost 

= tr[^Q (Q + gJrG^)] + V*(Gq) (5.4.12) 

The Principle of Optimality states that these two minimizations 

* T-1 

result in the minimizing sequence (^^).j._q 


for Problem AE. 



no 


Prom [Athans,411, "the Hamiltonian for the minimization over 
{l,2, . . - ,T-1> is 

^'^i,t^i=0' ' ^t ^ 

L 

= tr IT. . {Q + ) 

"^t-1 "" 


+ tr 


L3^ 


— > p . .IT. (A+B . G ^ ) E . (A+B . G ^ ) 1 S . ^ . 

TTj ^ 3-t -i,t ' D-t y-3,t+i 


for t e {l/2,3, . . . ^T-l} 

where the costate matrix is {S . . , ^ 

-],t+l 3=0 


(5.4.13) 


Remark ; We have now foirmulated Problem AE-1, which minimizes the accumu- 
lated cost over the interval {1,2 ,...,t} with respect to the sequence 
T“1 

1 using the matrix minimum principle and results in the optimum 
* 

cost, given G^ , V (^q )• Problem AE— 2 is then the minimization of 
equation (5.4.12) over • 
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5.5 The Necessary Conditions. 


The matrix minimum principle yields necessary conditions which 
an optimum must satisfy. There are two conditions of importance. 
{The third condition yields equation (5.4.6)). 

From the necessary condition for the costate, 

3h 


* 

^i,t 


. . 

-i,t 


the propogation of ^ . backward in time is derived . 

i , t 

I T 

^i,t = ^ 

IT • T T 

+ > , p . . - — [A S. ^,,A + G^B.S. 

^ IT. — -D,t+i- -t— : -:,t+i-:— t 


(5.5.1) 


+ a"^S . B . + G^ bTs . A] [. (5.5.2) 

— — 3,t+l— 3— t — t — 3— 3,t+l— j 

This equation is well-defined for any sequence t>0. 

The cost V of using this arbitrary sequence over the interval 

{1,2,...,t} is given by 






(5.-5.3) 


The total cost over the interval {0,1 ,...,t} using this sequence is 


= tr 
T 


S . , S . , 

— 1,1 — 1,1 


+ tr[(Q + G|^ RGq)Eq ] 


(5.5.4) 


= tr 


(WB . )’'S . j . I „ 




= tr 




(5.5.5) 


(5.5.6) 
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Define 




(5.5.7) 


Then from equations (5.5.6) and (5.5.7) 


J„ = tr[Z„S^] 
T — 0 — 0 


T-1 


Thus, the cost of a given sequence ( length T is 

(Gq ,G^ ] 

I 

For future reference , define the matrix ^ . by 

s’ A ^i,t 

— i,t TT. 

^t-1 


(5.5.8) 


(5.5.9) 


(5.5.10) 


and note that equation (5.5.2) becomes 


T T T T 

S_ = Q+G^RG^+ S.^^a + G^B.S.^^-B.G 

— i,t ^ -t 1 Z-J ^31 — — 3,t+l— — t — ] — 3,t+l— 3 “ 


+ a'^S . B . G^+ gJ^bTs . _ , A] (5.5.11) 

— — j,t+l— 3 — t — t— j— j,t+l — 


From the Hamiltonian minimization necessary condition 
3 h 


3£t 


= 0 


(5.5.12) 


the following relation between E . , S . , , and G is obtained. 

— l,t — 3,t+l — t 


^ 71. 2 . ^ 




3=0 j 


t •- 


[bTs . B . G^ + bTs . A] p..7r. Z . ^ 

— D/t+1— D— t — 3— D,t+1— ^^^ 32 . i i/t 

1— U t— 1 


(5.5.13) 
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Remark: At this point, a two-point boundary value problem has been 

defined with the constraint (5.5.13) relating equations (5.5,2) and 
(5.4.6). Equation (5.5.13) is not explicitly solvable for 
because ^ . cannot be factored out of the sum over j ; thus , it cannot 
be used as a substitution rule in the other two equations. At this 
time, the solution for appears intractable. Thus, although necessary 

4e 

conditions for the existence of , the minimizing gain, have been 

ic 

established, they ‘do not readily allow for the solution of , and 
certainly do not admit a closed-form expression. 



5.6 Problem B: The Non - Swi t ching Solution . 

Although the methodology presented in Section 4 yields the 
necessary conditions for an optimum, these conditions are not analyti- 
cally illuminating. In this section, a second optimization problem 
is formulated. An equivalent formulation was presented in [Birdwell & 
Athens, 40]. The solution will admit a closed form expression for u^ . 
Although this solution is not the optimal .solution for the first 
problem, in that this solution does not necessarily satisfy the neces- 
sary conditions for problem AE, it will be proved that the two solu- 
tions are equivalent in the sense that for the steady-state solutions, 
as defined in Section 7, either both solutions stabilize the system, 
or neither one stabilizes the system. Even better, it will be proved 
that the steady-state solutions to the problems are identical. 

For the system (5.3.1), the objective is to minimize at each time 

t the weighted sum, with respect to expected costs-to-go, 

* 

given the control (x.^) ^nd forT>t, and given 

that the structure at time t-1 was k(t-l) = i, for each i. 

Formally, let C be the expected cost-to-go, given 21 .f. ' — t ' 

k(t-l) at time t be defined, as 

C(x , u . , k(t-l) , t) = x^Qx.^ + u^R u^+ 

— t — t — t— — t — t— — t 

\rt) ,k(t) ,t+l) I k(t-l)l (5.6.1) 

* 


where * denotes the optimum value, and u ^ is computed as 


u = arg min 


(5.6.2) 


T 

= arg min TT £(t) 


(5.6.3) 
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and 


C (x ^ ,k(t-l) ,t) = C(x^ ,k{t“l) ,t) 


( 5 . 


where 


C(t) = 


C(x^ fk(t-l)=0,t) 




( 5 . 


and 


C(T) = C*(T) = 


T ^ 




( 5 , 


Thus, the problem is 


S T T 

ir. ^ 514 - 0 . 21 ,- ■*■ H. 

. - - 


+ ECC»(x ,k(t),t+l) k(t-l]=ill 
t+i 


( 5 , 


mxn / ^ TT . 


T T 


+ >,p..C*(Ax +B.u , j ,t+l) 

]i 1 — 1 — t 


< 5 . 


From the formulation, is non-learning in that it depends only 

TT^ ^ for its knowledge of the past. Let C* be of the form 


C*(x^,k(t-l),t) = 


( 5 . 


Then for t = T, 


^k,T = 2 


( 5 . 


6 . 4 ) 


6 . 5 ) 


6 . 6 ) 


6 . 7 ) 


6 . 8 ) 

on 


6 . 9 ) 


6 . 10 ) 
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And equation (5.6.8) becomes 


S T T 

IT. x^Qx + u^Ru. 

1. , _-t— — t — t 1 

t-1 

— t •*- 


p . . (A X ^ + B . u ^ ) ^ - (A 

— t -j —t —j ,t+i — 




(5.6.11) 


At the minimum, differentiating (5.6.11) with respect to , we 


obtain 


T * T * 

7T. RU.+ 7P..(B.S. B . U . t B . S . . . ^ A X , 

1 ^ ;]1' -t -D-D/t+1 1 


Solving for , 


(5.5.12) 




(5.6.13) 


and hence the gain matrix is given by 


^ . 1 -1 L 

R+> > tt.B.S. 

- ^ 3^-l-D,t+l-: J ^ 3^-3 -3, 




"k * 

where t 

From (5.6.11) and (5.6.4), 

y ' IT I * 

^t^k,t^t == 2 -^^t^^t 


t+1- 


(5.6.14) 




(5.6.15) 


or, since (5.6.15) holds for all / 

^ *ip * 

S. = Q + G.RG^ 

— k,t — — t 1 


T ’ T • * *T T • 

P -1, S . A + A S . B . G ^ + G ^ B . S . A 
•^]k — — j,t+l— — — 3,t+l— 3— t — t— 3— D,t+1 — 

+ G*B.s! .,,B.G^) 

— t - j —3 ,t+l ~3 — t 


(5.6.16) 
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Thus, (5.6.16) proves by induction that equation (5.6.9) is valid. 
Note that equations (5.6.16) and (5.5.11) are identical. 

•k 

Therefore, the unconditional cost of G , t=0,l, . . . ,T-1, is, from 


— t 


(5.5.9) 


Jt = 


(5.6.17) 


which in this case is simply 


St-i'S-o 


(5.6.18) 


The matrices are called the non-switching , or non-learning gains, 

* 

and will hereafter be denoted G . The label G ^ will be reserved for 

—ns. — t 

t 

the solution to equation (5.5.13). The optimal value of the cost-to-go 
at time t=0 for this problem will be called the non- switching cost index , 
and is given by 


us„ 


T 


= X 


S IT. x^S . . X. + x^ (Q + g"^ RG )x 

iQ-l-r,l-l -0 - -ns^— nsQ -0 

o / . . G )'^S . , (A+B . G ) 


+ Q + G RG 

- -ns^^~-nSQ 


^0 


(5.6.19) 


(5.6.20) 


Note that if G = G. for all time (i-e., if the solutions to the 
-nsQ -t 

optimal control gain problem and to the non-switching control problem 


are the same, then E [J ] = J 

X ns T 

T 
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Sximinary ; In this Section, the non-switching, or non-learning, gains 
have been derived. These gains are called non-switching or non-learning 
because they do not depend on the past trajectory of and u ^ , but 

only on the initial probability vector over 1/ 2 Lq • further 

shown that if the solutions to Problems AE and B were identical, then 

E IJ ] = J (5.6.21) 

X „ ns T 

— 0 T 
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5.7 Stability and the Steady-State Solutions . 

In this Section, the concept of stability for this class of 
systems will be precisely defined. From this, a natural concept of a 
steady-state solution to Problems 2\E and B will be given, and a very 
strong result relating the solutions to the two problems will be 
proved . 


5.7.1 Stability and Cost-Stability. 

For this class of systems , two definitions of stability will be 
tendered. The first is the usual definition of mean-square stability; 
the second definition, that of cost-stability, has a strong relation to 
the existence of solutions to the infinite time versions of Problems AE 
and B . 


Definition 1 : (Stability) . ^ is a constant stabilizing gain if and 

only if the resulting system given by equation (5.3.1) and repeated here 




~t+i 

is mean-square steible: 
E[jc^x^ ] -»■ 0 


as t 


(5.3.1) 


(5.7.1) 


Definition 2 : (Cost-Stability). The system (5.3.1) is cost-stable, 

if and only if the scalar random variable 


CO 

S 

^=n 




+ U 



< CO 


(5.7.2) 


with probability one. 
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5,7.2 Definition of the • Infinite-Time Cost. 

In this research, the infinite-time problem is defined as a 
minimization of 

J = lim (5.7.3) 

where is the cost function for the corresponding finite-time problem. 
The sequences which solve these infinite-time versions of Problems AE 

* oo oo 

and B are (G^ ^ and (G ^ t respectively, when a solution exists. 

— t t=(j — ^®t 

A Solution will exist if there exists a sequence of gains for which the 
limit in equation (5.7.3) exists. This definition of the infinite-time 
problem is chosen rather than the definition requiring a minimization 
of the average cost per unit time 

J = lim ij (5.7.4) 

T-«® 


because there is a direct correlation between the boundedness of J 

T 

over all T for a constant sequence of gains ^ and mean square stability 
of the system (5.3.1). It is necessary, however, to prove that the 
set of problems for which is bounded for some sequence of gains is 
not vacuous. This fact is demonstrated by any of the convergent non- 
switching gain examples in Section 10. 

As further demonstration of the validity of using equation (5.7.3) , 
note that if 0 < < <», then the cost per unit time has a non-zero 

steady-state value, which implies that the system (5.3.1) is not mean- 
square stable since 


= tr[E (Q + RG )] 

— ss — — ss ss 


(5.7.5) 


where I and G are the steady-state values of Z and G . , 
— ss — ss — t — t 

T 

they exist, and, since Q + G RG is positive definite, Z 

— — ss ss — ss 


when 
¥ 0 . 
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5.7.3 Bounded Cost and Mean-Square Stability. 

In choosing equation (5.7.3) as the basis for the definition of an 
infinite-time problem/ a major requirement was that the existence of 
an infinite-time solution, namely of a sequence of gains which yields a 
finite cost in equation (5.7.3), imply mean-square stability. For 
the case where the sequence is constant, the following result is 
proved. 

CO 

Theorem 1 ; A constant sequence of gains (G^)^_q is mean-square stabiliz- 
ing if and only if there exists a bound B < °° such that 

< B for all T (5.7.6) 

Proof : See Appendix 5.1. 

TO CO 

Remark: For a sequence (G ) <B<«»WT implies CG ) ^ is 

— t t— 0 t ' — t t=0 

CO 

mean-square stabilizing, but (^^) ^_^ mean-square stabilizing does not 
imply is bounded for all T. 

Proof ; See Appendix 5.2. 

5.7.4 Cost-Stability. 

As yet, the definition of cost-stability has not been utilized. 

In this Subsection, it will be shown that the system described by 

00 

equation (5.3.1) is cost-stabilized by a sequence of gains (G^^)^ Qif and 
only if J is finite-valued for this sequence. One direction of this 
result is proved in the following theorem. 

00 

Theorem 2 ; Any sequence (G^)^_Qfor which j < cost - stabil±zes (5.3,1) 
with probability 1. 

Proof; See Appendix 5.3. 
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The other direction of this result is obvious ; If a sequence 

CO 

(G^)^_Qis cost-stabilizing with probability one, then the random cost, 
given by equation (5.7.2), is finite except on a set of structural 
trajectories of measure zero. (The appropriate measure on this set is 
given in the proof to Theorem 2.) Since the expected cost J is the 
integral of equation (5.7.2) with respect to the probability measure 
on the set of structural trajectories (see Appendix 5.3) , then J is 
finite. 

Thus, the cost-stability and the existence of an infinite- time 
solution are equivalent. 

5.7.5 Equivalence of Problems AE and B. 

The first major result of this Chapter will now be stated. This 
result establishes a strong equivalence between the solutions to 
Problems AE and B. 

00 

Theorem 3: A cost-stabilizing solution (G ) „ exists if and only if 

—ns t=0 

. t 

* OO 

there exists a cost-stabilizing solution , assuming tt^> 0 for 

all i and Zq>0. 

Proof ; See Appendix 5.4, 

Remark 1 ; This result provides a computationally feasible methodo— 

03 

logy for arriving at a sequence of gains (G ) „ which cost-stabilize 

— ns^ t=0 

the original system (5-3.1) with probability 1, whenever such a se- 
quence exists. The coupled matrix equations of Problem B (5.6.16) can 
be iterated backward in time. If the weighted sura with respect to the 
ergodic distribution 7T converges, then the resulting sequence of gains 
cost-stabilises the system (5.3,1) with probability one. 



5.7.6 The Steady-State Solution. 
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A steady-state solution to optimization Problems AE and B can 
exist only if there exists a steady-state probability distribution tt 
over the set of possible configurations indexed by I such that 

TT = P TT (5.7.7) 

and 

lira TT = TT (5.7.8) 

— t — 

From equation (5.7.7), it is apparent that for TT^ to exist, the matrix 
P must have an eigenvalue at 1, and Ijl must be in the subspace spanned 
by the eigenvectors of IP corresponding to that eigenvalue. The fol- 
lowing lemma states precisely when ir exists. 


Ileimna 1 : ^ exists if and only if one of the following three conditions 

is satisfied for each diagonal element of the Jordan normal form ^ 
of P, where 


P = T A T 


-1 


(5.7.9) 


a 6 
0 o 


Oil 3i 
0i2 


A = 


L-1 




6^ = 0 or 1 


(5.7.10) 


For each i 
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i) |a. I < 1 

ii) a. = 1 
1 

ill) 1/ 7^ 1, (T ^ tTqI = 0 

Proof : Obvious . 


5. 7. 6.1 Steady-State Solution to Problem AE. 

•k k 

Note that for Problem AE, initially, the gains 

k k 

will depend on and near the final time, the gains ... G , 

U T 


will depend on a time-varying Thus, the steady-state solution for 

Problem AE is defined as the limiting solution to equations (5.4.6) 
(5.5.2) and (5.5.13) at time t, first as T^ and then as t->«>, if this 
limit exist. The steady-state values for r and , when 

they exist, satisfy the following equations: 


I. 

r] 


1_ 

IT . 
D 



(A+B . G ) 

D “ 


E. ■ (A+E.G)*^ 
“1 J~ 


(5. 7'. 11) 


S. =ir. 
—1 1 


T IT T T T 

Q + GRG+ > p..— [A S.A + G B.S.B.G + A S.B.G 

J- 3-J- 




T T 

+ G. B.S .A] 

3-3- 


(5.7.12) 


IT T 

0 = RG > IT. E. + > [B.S.B.G + B.S. A] > p . . IT . E. 

“ ^ ^ l-3^~3- -3-1- ^ 31 x-i 


(5.7.13) 


which are the limit of equations (5.4.6) , (5.5.2) , and (5.5.13), given 

k 

that the limiting solution E. and G exist, where tt satisfies 

— 3 ^ ^ ^ 

equations (5.7.7) and (5.7,8). The cost of this steady— state solution 
is 


J = lim J 

T 


(5.7.14) 


as in equation (5.7.3). 
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5. 7. 6 . 2 Steady-State Solution to Problem B. 

The solution to Problem B depends on its past only through the 
probability distribution Tr(t) over the structure index set I. 

Therefore / to develop the steady- state solution, let the initial pro- 
bability distribution TT ^ equal the steady-state value 7T from equations 
(5.7.7) and (5.7.8). Then the steady-state solution can be defined as 

the limit, when it exist, of the gain G .calculated for the problem 

— ns 

ending at time T, and of the solutions to the coupled Riccati-like 
equations (5.6.16), , as the final time approaches infinite. Let 

X / u 

I 

G (T) and S . „ (T) be the solutions at time zero for Problem B with 
— ns^ —1,0 

final time T. Then 


G 

—ns 


lim G (T) 
T-xx. 


t 1 

S_ . = lim S. (T) , i e I 
^ r , u 


(5.7.15) 

(5.7.16) 



T T ' 

+ G B..S.B.G ) 
— ns— j - 3 - 3 — ns 


(5.7.18) 


The cost of this steady-state solution, given x, is, when the limit 


exists 
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J 

ns 


lim J 

ns 

T 


T 


X 



t 

TT.S. X 
1—1 — 


(5.7.19) 


5.7.7 The Possibility of Limit Cycles. 

The discussions in the last Section do not rule out the possibi- 
lity of limit cycles in an infinite-time solution. In Problem B, 
the expected cost is directly computable from a set of coupled Riccati- 
like equations (5.6.16), as is the .non-switching gain (5.6.14). If' 
these coupled matrix equations converge whenever the solution is 
bounded, then the non-switching gain is always directly computable when 
it exists. Boundedness implies convergence of the expected cost 
(Lemma 2) ; however, the possibility of the existence of a limit cycle 
in the solution to equation (5.6.16) is not ruled out. It is con- 
jectured, but not proved, that such a limit cycle cannot exist. 

* 

Lemma 2 ; If the expected cost for Problem A is bounded, then it 
converges . 

Proof i See Appendix 5.5. 

* 

Since E [J ] = J , J 

x ns T ns 

T T 


also converges- 
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5.8 Equality of G and G 
ns 


In this Section it will be shown that when a steady-state G and 

— ns 

* * * 

G exist, with finite cost J and J , the gains are equal. This 
— ns 


result is extremely important in that it yields a method of calculating 
the steady-state solution to a two-point boundary value problem as the 
limiting solution to an equivalent (in the steady-state) single boundary 
value problooi. it is taken as a working hypothesis in this Saction that 
both problems have a steady-state solution and that the ergodic distribu- 
tions of and , for all i, exist. Then the steady-state cost of the 
optimal problem is 


J*. = tr[^^ (Q + G*RG* )] + 


ss 


-O-.-O 


* * 

tr[S . S . 
— 1—1 


For any constant gain ^ for which the limits exist, 
be 


) (5.8.1) 

the value would 


Jss^-^ "" tr[^Q(Q + G RG)] + ^ tr (G) S ^ (G^) ] (5.8.2) 


= tr[^^ (£ + £ 


RG )] + tr[ (A+B . G)2 (A+B . G)’^S . (G) ] 

^ ~- i -_0 -1 1 


= tr 


2„Io + G^RG + > (A+B . G)^S . (G) (A+B . G) 

”°r ^ "" 


(5.8.3) 


(5.8.4) 


Similarly, equation (5.8.1) becomes 


J = tr 
ss 


i **]! * * rp * * ) 

2 -*° 0^^0 * g i°o > i = 0 >} 


(5.8.5) 


For the non- switching, or non-learning problem, the steady-state cost 
for any ^ for which the converge is, given x.q 
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J (G) = X^(Q+G^RG)x^+E 

ns — —0 — — 0 —0 

ss 


(5.8.6) 


= iio ^2 


+ TT. (A+B . G)^S . (G) (A+B . G)x„ 

— 0 1 X— —X — X 0 


(5.8. 7) 


Taking ej^ectations with respect to x ^ ^ 


ns — 
L ss ■ 




+ > tr [Z - (A+B . G) S . (G) (A+B , G) ] 

^ „0 _„1„ _ __i_ 


(5.8.8) 


E (G) = J^(G) 

L "^s-J 


(5.8. 9) 


Thus, the costs are equivalent for any ^ for which the equations 
converge . 


By Lemma 3, if the non-switching expected cost is bounded for a single 
G, then the equations converge; i.e., there can be no limit cycle. 


Lemma 3 ; For a given gain ^ , if the expected cost (^) is bounded 
then it converges. 

Proof: See Appendix 5.6. 


Thus, either equation (5.8.9) holds, or both costs are infinite. There- 
fore, if the cost is finite for any single then there exists a 
which minimizes both costs. Furthermore, given that (T) converges, 

G (T) + G . as T +- CO. This result with an extension is stated in 


Theorem 4. 



129 


* I 

Theorem 4: Assume the values G . (T) , G (T) , s . (T) , S . (T) , and 

— t — ns^ — lyt — i,h 


^ converge. Then 


A) G (T) -> G ^ as T which minimizes equation (5.8.9J. 

— ns^ —opt ^ 

•k 

B) G = G , where G is the steady-state value of G (T) , 

— ns — —ns 


* . * 
and G is the steady-state value of G (T) : 

— — t 

lim lira G (T) = G (5.8,10) 

■t-X» rp-X)0 

Proof : See Appendix 5.7. 


Discussion : The result of Theorem 4 B) gives a direct computational 

* 

procedure for calculating the, optimal steady- state gain ^ as the 


limiting gain • There are, however, still some open questions 
concerning the existence of limit cycles in the calculation of G^^ 

OO 

Theorem 3, however, guarantees cost-stahility using (G‘ ). „ if a 

—ns^ t=0 


cost~stabilizing sequence of gains exists. 
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5 . 9 Robustness . 

The original problem (Problem A) can be formulated in such a way 

00 

that the sequence will cost— stabilize a set, of linear systems 

n s ^ "C u 

with different actuator structures individually whenever such a stabiliz- 
ing or robust gain exists. 

Definition 3 ; A gain ^ is robust if 

is stable for all k. This is the same as requiring the matrix (A+B 

iC 

to have eigenvalues inside the unit circle for all k. 


Corollary 1; For the set of L+1 systems 


X . - = Ax^ + B, u^ 

— 1+1 1 — k— t 


(5.9.2) 


with 


P = I 


(5.9.3) 


7T. 

3 


1 

L+1 


(5.9.4) 


00 


if a robust gain exists, then (^^^^ ^t— 0 ^ stabilizing sequence for 

(5.9.1) for each k, and if the gains ^ (T) converge, then G is a 


robust gain. 

Proof ; For the expected cost to be finite , for any ^ ^ must be 
robust, since each structure is equally likely and no structural changes 

ie OQ 

can occur. Therefore, if a robust G exists, then certainly (G^)^ ^ 

— — t t=0 

. 00 
will be stabilizing, and by Theorem 3, so will (G Also, if 

— nSj^ t=0 


G (T) converges as T the G will be robust since it will have 

^ —ns 
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finite cost J (G ), which implies stability, in this case, for all 
k e I. 

Q.E.D. 

Discussion ; With Corollary 1, a specific existence problem for robust 
linear gains is solved. Existence of a robust gain is made equivalent 
to the existence of a finite cost infinite-time solution to Problem B, 
which is readily computable from equations (5.6.14) and (5.6.16) . 



5 . 10 Examples. 


In this section, two examples are presented to illustrate the 
non-switching gain computational methodology. Example 5.1 is ana- 
logous to Example 3.1 of Chapter 3; it demonstrates the effect of 
component reliability on system stabilizability with a non-switch- 
ing- gain control law. The first case of Example 5.1 is not conver- 
gent; -the second case is convergent. The only difference between 
the two cases is the reliability of the actuators- Case i) corresponds 
to Case ii) of Example 3.1; Case ii) corresponds to Case iii) of 
Example 3.1. Neither case results in a robust control law, but ro- 
bustness is not possible because the system is uncontrollable in 
structural state 3. As an aside, it is interesting that the "optimal" 
non-switching gain in Case i) ignores state x^; the system is decoupled 
in that there is no interaction between and x^ • Since state 
has stable dynamics, and the dynamics of state x^ are unstable, the 
entire control effect is concentrated on state x^- 

The computer routines which are used in the calculation of the 
non-switching gain solution are listed in the Appendix. The primary 
subroutine is AIM; it calls WEIGHT. Any other routines which are 
used are from the standard ESL subroutine library. 



133 


Example 5.1; 


A = 


^0 = 


= 


.71828 

0.0 




.0 

.36788 

_ 




1.71828 

1.71828 

^1 = 

0.0 

1.71828' 

-.63212 

.63212 


0.0 

.63212 

[l. 71828 

o.ol 


o 

• 

o 

o.ol 


••63212 0.0 


^3 


0.0 0.0 



14. 

8." 


'l.O 

o.o' 

= 



R =: 




8. 

6. 


0.0 

1.0 


p = 


1-2p^-HP^ 


Pf-Pf 


Pf-Pf 


P. 


(l-Pf)P^ 


l-Pjr“P +P^P 
f r f r 


PrPf 


(l-P^)Pf 


(l-Pf)Pr 


PrPf 


1-P--P +P^P 


Cl-Pr>Pf 


P,(l-P,) 


P„{1“P„) 
r r 


l-2p +p 
r r 


The system is 

= Ax^ + B 


kCt) — t 


-t+1 1 

k(t) e {o,l,2,3} 

The cost to be minimized is 


2-t ■ '=l,t *2,P 


J = E 


Eit2»t^4£S,tl£ 


t=o 
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Example 5,1, Case i) 

- - 



p^ — Pj. — 

.81 




.09 


""l 

IT = 


= 



.09 


. '"2 


H 

O 




L .. 


3 


Non-Co nvergent; but gain converges at 


G — 
•“ns 


-1.246 

-1.039 


0.0 


0.0 


J 


Stability; 


Configuration Stable 



no 

1 ) 

yes 

2 (B,) 

yes 

3 ( 83 ) 

no 


Interpretation: The coupled Riccati equations are unbounded. Note 

that since state has stable dynamics, the convergent non-switching 
gain concentrates on stabilizing x^, which is open-loop unstable. 

From the above stability table , the control law 


stabilizes only configuration states 1 and 2 ; since the configuration 
has a high probability of being in state 0 (unstable), the cost diverges. 



.98 
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Example 5.1, Case ii) 


= - 01 , 




.9799 


■^0 

.009999 


^1 

.009999 


^2 

.0001020 


^3 





Convergent Coupled Riccati Equations. 


G 

~ns 


-.7563 .1266 

-.8070 -.1784 


Stability; 


Configuration Stable 

0 (B Q ) yes 

1 ) no 

2 (B 2 ) no 

3 (IB 2 ) no 


Interpretation; With more reliable actuators, the non- switching gain 
expends less force on the stabilization of configuration states 1 and 2 
(unstable) ; since configuration state 0 is stabilized, and the system 
has a (relatively) higher probability of being in configuration state 0 
than in Case i) , the non-switching coupled Riccati equations converge, 
resulting, in a finite cost. 
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Example 5.2 uses the same system dynamics as in Example 5.1; 
however, only structures 0,1 and 2 (the controllable structures) are 
considered. The configuration dynamics are modeled as being in any 
structural state with equal probability of occurance initially and 
remaining in that state forever; this model is illustrated graphically 
in Figure 5.1. 

The state dynamics are 


it+l 


K(t) e {0,1,2} 

The cost to be minimized is 


^t = ^^l,t ^2,t^ 


J = E 


t=0 


T T I 

X . Qx . + U . Ru . 7T 
t t t t 


The non-switching methodology yields a robust control law of the 
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76265AW030 



Figure 5.1: -Markov transition probabilities for Example 5.2. 



Example 5.2; 


2.71828 


A = 


O.-O 


0.0 

.3679 



1.71828 

1.71828’ 


0-0 


1-71828' 

il 

o 

-.63212 

.63212 

^1 = 

0.0 


.63212 



“1 









1.71828 

o 

» 

o 



1. 

0. 

0. 


-2 

-.63212 

0.0 


P = 

0. 

1. 

0. 







0. 

0. 

1. 



Convergent : 


G 

— ns 


-1.089 

-1.028 


-.008413 

-.01444 



112.8 

8.992 


8.992 

6.835 



Stability: 


Conf igura t ion 


Stable 


0 (B„) 

yes 

1 (B^ ) 

yes 

2 (B^) 

yes 


Robust : yes 
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Riccati Solution: 

, [l09.8 9.030 

S = 

° 9.030 6.821 



114.3 

6.285 


6.285 

6.836 


^2 


114.4 

11.66 


11.66 

6.849 



14D 

The non-switching solution converges for the system in Example 5.2, 
and the three resulting configurations are stabilized. Therefore 
is a robust gain. Had the solution not converged, by Corollary 1 of 
Section 9, no robust gain would exist. 

The apriori expected cost (before the configuration state is 
known) is, given x 
T 

J = X C X 
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5 . 11 Sunimary . 

In this Chapter, an optimization problem was defined on linear 
systems with variable actuator configurations and quadratic cost criteria 
The objective of this approach was to compute apriori a sequence of 
gains to be used in linear feedback control which do not depend on 
any on-line information about the process. These gains were to 
both stabilize the overall system, accounting for the various possible 
structures and minimize the expected value of the quadratic cost crite- 
rion, where the expectation is taken over the possible sequences of 
actuator configurations. This solution depends on both the perfor- 
mance, and on the reliability of the various structures, as represented 
by the Markov transition probabilities between structures. 

The matrix minimum principle [Athens, 41] was used to establish the 
necessary conditions for optimality of a solution to an equivalent 
deterministic problem to that described above, known as Problem AE in 
the Chapter. These conditions unfortunately do not yield an analytic 
solution for the gain sequence, but instead yielded an ill-posed two- 
point boundary value problem which must be solved numerically (Section 5) 
Therefore, a second problem (Problem B) was formulated which was solvable 
analytically using dynamic programming (Section 6) . This solution has 
identical cost-stabilizing properties to the solution of Problem AE, 
but has the advantage of being directly computable. 

The steady-state solutions to the infinite-time versions of both 
problems were defined, when they exist, and it was proved that, in addi- 
tion to the equivalent stabilizing property of the two solutions, the 
steady-state values are identical, and this value is the same as the 
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optimal constant gain which minimizes the expected cost over the infinite 
time interval . 

In addition, the general robustness question of when one gain can 
stabilize a set of linear systems with different actuator configurations 
was formulated in the context of Problem A and was solved by Problem B. 
Thus, a test for when a robust gain exists can be performed by iterating 
a set of coupled matrix Riccati-like equations and testing for converg- 
ence of a function of the solutions. If, in addition, the individual 
solutions converge, then the robust gain which minimizes the expected 
quadratic cost index can be calculated directly. It was noted that the 
extension to systems with variable dynamics (variations in A) , as well 
as variable actuator stmcture, is trivial as long as the dimension of 
the state is constant. 

The major applications of this work are in the calculation of a 
robust gain for a set of linear systems and in the calculations of 
stabilizing gains for systems with variable structure, such as occurs in 
failure, repair, or reconfiguration. A second application will be 
covered in the next Chapter and involves using these calculations in a 
computer-aided design procedure for the determination of the relative 
effectiveness of various redundant component configurations. 
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CHAPTER 6 

COMPUTER-AIDED DESIGN 


6.1 Introduction . 

In this Chapter/ two specific applications of the non-switching 
gain methodology to computer-aided design are presented. Example 6.1 
illustrates the usefulness of the non-switching gain methodology in 
the selection of an actuator design. Five possible designs are 
analyzed using the non-switching gain calculations as a basis for ranking 
the designs with respect to their expected performance. Example 6.2 
compares two actuators, of which one is more reliable, but less 
effective (in that it incurs a greater cost for the same action) than 
the other. Three cases with various actuator reliabilities are presented 
as a study of the trade-off between actuator reliability and effective- 
ness. 

These two examples are intended to demonstrate the usefulness of 
the non-switching gain methodology in design studies. No general method- 
ology for computer-aided design using the results presented in this 
report is presented. Instead, tools are presented which can be used in 
the computer-aided design of system configurations. 

6 . 2 The Design Decision . 

A designer often has many means of achieving a desired goal; 
however, no unified methodology exists which can be used to choose a 
given design that is "better" than any other. At best, a set of tools 
can be developed which are applicable to specific situations and classes 
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of systems. Of these tools, all that are presently available evaluate 


a system either on the basis of performance or on the basis of reliabil- 
ity. The methodologies described in this report optimize a performance 
index which depends on both system reliability and system performance. 
Therefore, it is logical to apply these methodologies to the computer- 
aided design of system configurations. 

Example 6.1 is an aid in the design of a linear system for which the 
state dynamics are fixed, but the actuator configuration is to be at 
most two actuators (one level of either component or functional redundancy) 
chosen from two types of actuators. The system in Example 6.1 is de- 
fined by 


-t+1 


Ax .+ B, . u^ 
1 — k (t) — t 


( 6 . 2 . 1 > 


k(t)ei ■ (6.2.2) 

where x 4 - “ 4 -'- Cases i) and ii) , I = {o,l} ; 

in Cases iii) , iv) , v) , I = {o,l,2,3}. The cost to be minimized is 

x"^Qx^+ u^Ru^ I IT (6.2.3) 



The cost of each actuator (labeled b ^ and b ^ ) is to~be the quadratic 
cost incurred by the control input to that actuator. These costs are 
represented by the quadratic weights r^ and r^, respectively, and are 
equal in Example 6.1. The actuators act on different states of the 


system; actuator b^ applies the control force to state x^, while 
applies the control force to state x^- Each actuator can fail to an 
actuator with zero gain, Repair constitutes replacement of the 
failed component with a new actuator, identical to the original ac- 
tuator. The repair action is modeled using a Markov transition pro- 
bability p^, the probability of repair per unit of time. The actuators 
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have identical probabilities of failure and repair per unit time, p^ and 
p^, respectively. The five possible actuator configurations are, in the 
order in which they are presented in Example 6.1, 


b" = 

l^ol 


(6.2.4) 




(6.2.5) 

= 

|l^0 i 

^ol 

(6.2.6) 

II 

1^1 1 


(6.2.7) 

b" = 


^ll 

(6.2.8) 


Configurations ^ and B have two-state configuration dynamics directly 

defined by the failure and repair probabilities per unit time. Con- 
3 4 5 

figurations B , ^ and B have four-state configuration dynamics re- 
presented graphically by Figure 3.2 of Chapter 3, Section 5. It is 
not immediately obvious from the configurations and the state dynamics 
which configuration is optimal. When a non-switching gain control is 
used, the expected steady-state cost, given by equation (5.7.3), is 
a measure of the expected performance of each configuration, and can be 
used to rank the five configurations in order of system effectiveness ■ 
System effectiveness is a measure of the expected performance of a 
system, taking into account all postulated modes of operation. There- 
fore, in Example 6.1, the non-switching gain and expected cost is com- 


puted for each of the five design conf iguraitons . 
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Example 6.1: 


2.0000 

.5000 

.5000 

A = 0.0 

0.0 

1.000 

0.0 

-1.000 

0.0 


1.0 0.0 0.0 

Q = 0.0 1.0 0.0 

0.0 0.0 1.0 



0 



Example 6;l case i) 


^0 = 

[b J = conf. 

, 0- 

(conf. 

^1 = 

0 = conf . 

. 1 

R = 


l-Pjr P 

f r 


.99 .98 

F — 

p, l~p 

f r 


.01 .02 


(conf. = configuration) 



[l.O 

0 

1 


IT - 


',9899 



.01010 


7i. 


Convergent Coupled Riccati Equations: 

■ = h„ 


■ ns 


,863 -.2582 -1.733 




182.5 

37.06 

57.93 


^0 

= 

37.06 

9.943 

12.32 




57.93 

12.32 

22.81 




188.6 

37.39 

60.09 


t 

-1 

= 

37.39 

9.961 

• 

12.44 




60.09 

12.44 

23.58 





182.6 

37.07 

57.95 


t 

ir.s . 

1— X 

= 

37.07 

9.943 

12.33 




57.95 

12.33 

22.82 


^ c 


Expected cost = x C^ x 
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Stability ; 


Configuration 


■0.(B„ )- 


KBi) 


Stable 

yes 

no 


Interpretation: The steady-state non-switching gain exists; it 

stabilizes configuration 0 (B.q) j ii^it does not stabilize configuration 
1 ' Since the probability of being in configuration 0 (stable) 

( TT^ } is much greater than the probability of being in configuration 1 
(unstable) (ir ^) , the system configuration is stabilized using the 
non-switching gain in the control law 


u ^ = G 
— t — ns— t 
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Example 

6,1 

Case ii) 

® n = 


= conf . 0 

— 0 

— 1 


^1 = 

0 

— conf . 1 


P = 


P 

ITf 

P* 1~P 

r 


R = r 


.98 

.02 

J 


.9899 


.01010 TT 


Convergent Coupled Riccati Eguations: 


G = -12.59 -1.484 -4.097 

— ns 


1035. 125.0 271.4 

125.0 18.84 33.04 

271.4 33.04 73.80 

1069. 129.0 282.6 

129.0 19.31 34.34 

282.6 34.34 77.43 

ri035. 125.0 271.6 


TT S . = 125.0 18.85 33.05 = C 

1—1 — 

271.6 33.05 73.83 


Expected cost = x C x 



stability: 


Configuration Stable 
-0 yes 

1 no 

Interpretation: The steady-state non-switching gain exists; it 

stabilizes configuration 0 ? but does not stabilize configuration 

1 Since the probability of being in configuration 0 {stable) 

( 7T^) is much greater than the probability of being in configuration 1 
(unstable) (7r^) , the system configuration is stabilized using the 
non-switching gain 


in the control law 



Example 6.1 Case iii) 


2-0 = 

~0 1 

^0 

- conf. 

0 

2-2 = ■ 

^ol 

£ 

= conf 


0 ! 

^0 

= conf. 

1 

2-3 " 

0 i 

£ 

= conf 


F'o °-°l 


o 

. • 

o 

0 

1 - 

1 

o 

u 

o 

• 

o 

— 1 


0.0 l.oj 


p = 


l-2p^+p^ 


Pj.(l-Pf) 

2 

Pf (1-Pf ) 

i-p^-pj+p^£ 

P^f 

p^li-p^) 

(1-p^) 


1-P^-Pj+P^j 

Pr'l-Pr' 

■ 2 
Pf 


I’£<^-Pr> 

2 

l-2p -fp 
r r 


I 

t 


.9801 

.9702 

.9702 

‘.9604 

i 

.0099 

.0198 

.0198 

f 

.0196 

.0099 

.0098 

.0098 

-.0196 

.0001 

.0002 

- .0002 

'. 0004 


r *1 

.9799 



.009999 


^2' 

.009999 



.0001020 


4 

■^4 - 



y 



Convergent Coupled Riccati Equations : 



'-2. 469 

-.1279 

- . 8983' 

G' 

— ns 




.-2.469 

-.1279 

-.8983, 






153.1 

32.81 

48.01 ■ 

■ 

^0 = 

32.81 

9.050 

10.92 


48.01 

10.92 

19.03 


- 




154.4 

32.88 

1 

CO 

CO 

f 

32.88 

9.054 

10.95 


48.48 

10.95 

19.20 





. 



-| 


154.4 

32.88 

48.48 

t 




^2 = 

32.88. 

9.054 

10.95 


'48.48 

10.95 

'19.20 


- 








155.8 

32.95 

48.96 

1 

^3 = 

32.95 

9.058 

.10.97 


48.96 

10.97 

19.38 



153.2 

32.82 

48.02 

32.82 

9.050 

< 

10.92 

48.02 

10.92 

19.04 


T 

Expected cost = x C x 
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Stability: 


Configuration Stable 

0 (B Q ) yes 

1 (B ) no 

2 (B 2 ) no 

3 ^ 


Interpretation: The steady-state non-switching gain exists; it 

stabilizes configuration 0 ) / i^ut does not stabilize configurations 

1,2, or 3. Since the probability of being in configuration 0 (stable) 
(tTq) is much greater than the probability of being in any other con- 
figuration (tt^, i=l,2 or 3) (unstable), the system configuration is 

stabilized using the non-swithcing gain G 

“ns 


u . = G X 
— t —ns — t 


in the control law 
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Example 6.1 Case iv) 


^0 

= 

^1 


= conf. 0 

^2 = 

^ll 

0_j 

= conf 

-1 

= 

0 1 

^1 

. = conf. 1 

^3 = 

0 i 

0 

= conf 

R = 

^1 

0.0 

o . o ' 
^1 - 

"l.O 

- 0,0 

o . o ' 

1,0 





P and IT are the same as for Case iii) . 
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Convergent Coupled Riccati Equations ; 




-6.097 -.7347 

-2.011 


G 

—ns 


-6.097 -.7347 

-2.011 

. 



762.2 

95.14 

195.1 


-0 

- 

95.14 

15.18 

24.64 




195.1 

24.64 

52.13 




768.7 

95.92 

197.3 


1 

^1 

= 

95.92 

15.27 

24.89 




197.3 

24.89 

52.83 




768.7 

95.92 

197.3 


1 

-2 

= 

95.92 

15.27 

24.89 




197.3 

24.89 

52.83 




775.3 

96.71 

199.5 




96.71 

15.36 

25.16 




199.5 

25.16 

53.55 





762.3 

95.15 

195.2 

s 

1 

TT.S . 
1—1 

= 

95.15 

15.18 

24.64 




195.2 

24.64 

52.14 


= C 


T 

Expected cost = 21 C x 
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Stability: 


Configuration 


Stable 


0 ) yes 

1 ) no 

2 (B ^2 ^ no 

3 ) no 


Interpretation: The steady-state non-switching gain exists; it- 

stabilizes configuration 0 <^q), -but does not stabilize configurations 
1, 2, or 3. Since the probability of being in configuration 0 (stable) 
(TT^) is much greater than the probability of being in any other con- 
figuration (tt^, i=l,2 or 3) (unstable), the system configuration is 
stabilized using the non-switching gain ^ 


-t 


G X 
—ns — 


t 


in the control law 
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Example 6.1 Case v) 



P and IT are the same as for Case iii) . 



Convergent Coupled Riccati Equations: 


G 

—ns 




1 





■“3-815 

-'.1312 

-1.106' 

.-2.956 

-.5815 

-1.486. 

126.5 

24.86 

32.32 

24.86 

7.066 

6.842- 

32.32 

6.842 

10.69 

ri28.4 

24.93 

32.88 

24.93 

7.096 

6.863 

00 

00 

CN 

ro 

6-863 

10.85 

127.3 

25.01 

32.72 

25.01 

7.097 

6.921 

32.72 

6.921 

10.89 


129.2 

25-08 

33-28 


25.08 

7.100 

6.942 


33.28 

■ 6.942 

11.05 



126.5 

24.86 

32.33 

= 

24-86 

7.067 

6.843 


32.33 

6.843 

10.69 


T 

Expected cost = x 


C x 
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Stability : 


Configuration Stable 


0 (B Q ) 


1 (Bi) 

2 (B , ) 


3 

Interpretation: The steady-state non-switching gain exists; it stabil- 
izes configuration 0 and 2 and ^2 ^ the probabilities of 

being in configuration 1 and 3 (B^and ) are small (7T^ and ir^) 
(unstable), the system configuration is stabilized during the non-switch- 
ing gain the control law 

= G X. 

— t — ns — t 
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From the results in Example 6.1, the design configurations are 
ranked as follows, where > is defined as "is better than". 


5 3 14 2 

B >B >B >B >B 


(6,2.9) 


One configuration is more desirable than another > B ) if 




> 0 


(negative definite) 


( 6 . 2 . 10 ) 


This criterion is reasonable; if >■ B^ r then the expected cost using 

*1 3 ^ 

design configuration is always less than that using ^ . If the left 
hand side of equation (6.2.10) is not negative definite, but is only 
semi-definite, then some other criterion must be used in addition to 
(6.2.10) to rank the various designs. For example, if one assumes a 
uniform distribution of the initial system state in the unit sphere, 

and if the elements of the diagonal of the left hand side of equation 
(6.2.10) are all non-positive, then the trace operator may be used as a 
ranking function. If the trace of the left hand side of equation (6.2.10) 
is negative, then ^ • If the left hand side of equation (6.2.10) is 

not semi-definite, then the designer must choose which of the state 
variables are most important in an effort to eliminate the ambiguity of 
equation (6.2.10). In Example 6.1, equation (6.2.10) alone is sufficient 
to rank the designs. 

The results stated in (6.2.9) are somewhat surprising. First, 

consider and ^ control input at time t using enters the 

T 

system dynamics in state x_, where x^ = [x. ^ x., x_ ^] . At time t+1, 

3 — t l,t 2,t 3,t 

the same control is applied to state x^ with a gain of .5? also, 

X = X . At time t+2, that control is again applied to state x^ 

with a gain of .5 . Now, consider 
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the same situation, but with b instead of b^ . In this case, at time 
t+1, the control is applied to state with a gain of .5, but 

= -x„ . . Therefore, at time t+2, the negative value of the original 

j) / t"p d. 2 / tl 

control is applied to state x^, thus partially cancelling the effect of 
the original input. The same process occurs using , but is delayed 
one time step; thus, the control affects state x^ positively one additional 
time step when is used. Because of the added effectiveness of b^ 

12.14 

over ^ 2 ^ / B ^ B^ / and in fact, B > B. Thus, even after accounting 

1 

for component reliability, configuration ^ , which has no component 

2 4 

redundancy is more desirable than configuration B or jB even though 
4 

configuration ^ employs one level of component redundancy. 

Using this reasoning, one would expect B^^ to be the optimal design 

choice; however, the example demonstrates that this is not the case. 

From G for Case iv) , note that the control which is applied to b _ 

—ns — O 

depends mostly on the unstable state x^, while more emphasis is given 
to states x^ and in the calculation of the control for actuator • 
Thus, actuator b^ acts partially to stabilize the dynamics of state x^, 
while actuator b^ acts partially to counteract the negative effects of 
the subsystem of states x^ and x^. This type of control action is an 
example of the use of functional redundancy, and is not possible with 
design configurations ^ or 

The non- switching gain analysis of the proposed design configura- 
tions yields information not only about the effect of various actuator 

configurations but also about the effect of component reliability on 

4 2 3 

the expected performance. Thus, B^ is more effective than B^ , and 

14 3 

is more effective than ^ ; B_ and ^ are versions of the configurations 



2 1 

^ and ^ , respectively, with one level of component redundancy.. Con- 
figuration ^ is an example of functional redundancy; both actuators 
provide control input to the same system, but are not identical components. 
Thus, the additional reliability of component redundancy contributes 
to ranking (6.2.9). The trade-off between system performance and system 


reliability will be further demonstrated in Section 3. 
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6.3 A Trade-Off of System Performance Versus Reliability . 

The non-switching gain methodology can be used to study the 

relative effects of actuator reliability and actuator effectiveness 

on expected system performance. If a designer has a choice between 

using a high reliability actuator rather than one with relatively low 

reliability, but with a higher effectiveness, on what basis can a 

decision be made? In Example 6.2, two actuators are considered. Each 

actuator may fail to an actuator of gain zero (0) and be repaired 

(replaced) . The probabilities of failure and repair are p^ and , 

i i 

where i=0 or 1 and refers to the actuator (b ^ ' respectively) . 

One actuator (b^ ) has good reliability, but the actuator gain is unity- 
A second ■ actuator (b^ ) has an actuator gain of ten (higher effective- 
ness) , and a lower reliability. If the actuators had the same relia- 
bility, then actuator b ^ would be preferable — it incurs a smaller cost 
for the same effect. In Case i) of Example 6.2, this reasoning is 
demonstrated numerically; the steady-state non-switching gain favors 
actuator b ^ (the second column of ) • (The two rows of the gain 
matrix are compared; the top row corresponds to actuator b^.) 

In Cases ii) and iii) of Example 6.2, the reliability of actuator 
b^ is lower than the reliability of actuator • in Case ii) the 
probability of failure per unit time of actuator b^ is five times 
greater than the probability of failiure per unit time of actuator b ^ ; 
in Case iii) , it is ten times greater. The probabilities of repair per 
unit time for actuator are also lower than for actuator . 
Therefore, actuator is significantly less reliable than actuator b^ 
Note that in Case ii) , the optimal non-switching steady-state controller 
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favors actuator by a gain factor of 2,5 - 2.6; in Case i) , actuator 
b^ is favored by a gain factor of 2.3. In Case iii) , actuator b^ is 
favored by a gain factor of 5.1. Tnus, the non-switching gain calcula- 
tions -can be- quite' sensitive to changes in component reliability. 

Although the configuration states are identical for all three Cases of 
Example 6.2, the configuration dynamics are modified by the changes in 
actuator reliability. The effect of modifications in actuator reliability 
on the non-switching steady-state gain and cost is pronounced. The 
s teady- s tat e gain is very sensitive to the actuator reliabilities ; the 
expected steady-state cost increases as the reliability decreases . A 
second effect demonstrated by Example 6.2 is interesting. In Case i), 
configuration state 2 is not stabilized by the non-switching gain. As 
the reliability of actuator decreases, the average steady-state 
probability that the configuration is state 2 (actuator failed, 
actuator operational) increases. Therefore, the non— switching gain 
solution must concentrate more effort on stabilizing configuration state 
2. Note that in Cases ii) and iii) , configuration state 2 is stabilized 
by the non-switching gain solution. It is interesting to note also that 
the non-switching gains in Cases ii) and iii) are robust with respect to 
configuration states 0, 1 and 2- (Configuration state 3 is uncontrolla- 
ble.) 


The system dynamics in Example 6.2 are 


— t-H 


= A X ^ 


+ 


— k(t) — t 


k(t) £ I 


(6.3.1) 

(6.3.2) 


where 


I = {o,l,2,3} and x = [x . x 




T 


The set { B . „ 

I — i -'1=0 


of configuration states is given in Example 6.2, The cost to be 



minimized is 
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( 6 . 3 . 3 ) 
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Example 6.2; 


A = 


2.0000 .5000 .5000 
0.0 0.0 1.0000 
0.0 - 1.0000 0.0 


^0 




Q = 


0-0 0.0 

0.0 0.0 

1.0 10.0 

0.0 0.0 


= conf. 0 


1 0.0 

0.0 

= 

1 1.0 

0.0 


1.0 

0-0 

0-0 

0,0 

1.0 

o 

• 

o 

0.0 

0.0 

1.0 


= conf. 2 B 


R = 


o 

o 

0. 

o 

o 

0. 

o 

d 

1 

10. 

o 

• 

o 

1 

0. 

O 

• 

o 

0. 

1 

o 

o 

0. 

5.0 

o 

o 

0-0 

5.0 


= conf . 1 


= conf. 3 


P = 


1-Pf -Pf +Pf Pf 
1 2 1 2 


p (1-p ) 

1 2 


p (1-p ) 

2 1 


^f 

1 2 


(1-Pf )P^ 

2 1 


1-P^ -Pf +P^ Pf 
1 2 1 2 


Pf P^ 

2 1 


(1-P^ )Pf 
1 2 


(1-Pf )P„ 

1 2 




'1 2 


P P 
r r 
1 2 


(1-p )p 
r r 
2 1 


1-P^r -P„ +Pf P^ (1-P^ >Pt- 
1 2 1 2 1 2 


(1-P„ )Pf 1-P„-P„+P„P„ 

2 1 1 2 1 2 
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P. = .01 

0 

= .01 p = .98 

f r 

1 1 


.9799 

.009999 

.009999 

.0001020 


Convergent: 



’-.2059 

-.01076 

-.07574 

G 

—ns 

-.4829 

-.02505 

-.1789 _ 


134.5 

30.06 

41.49 

30.06 

8.459 

9.981 

41.49 

9.981 

16.44 


134.5 

30.06 

41.49 

30.06 

8.459 

9.981 

41.49 

9.981 

16.44 


138.5 

30.27 

42.96 

30.27 

8.470 

10.06 

42.96 

10.06 

16.98 


138.5 30.27 42.97 


= 30.27 8.470 10.06 

42.97 10.06 16.98 


IT. S 
X — 1 


134.5 30.06 41.51 


30.06 8.459 9.9821 = c 


41.51 9.982 16.45 
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T 

Expected cost - x ^ 5S. 

Stability: 

Configuration Stable 
0 ) yes 

^ 1 ) y®® 

2 (B 2 ) . no 

3 (B^2 ) 

Interpretations : The system x.,_ =[A+B.G ]x is stable only for 

— t+i — — L — ns — t 

i=0 and 1. The probabilities of the configuration being in states 2 and 3 

(7T2 and tt^) are small; the system configuration is stabilized using the 

control gain G in the control law 
—ns 
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Pf - 

.01 p 

= .98 



”.9378 


^0 

0 

0 




.009212 




.05 p 

^1 

II 

o 

TT 


.05206 


^2 






.0005316 


'n'o 






- 


L_ 3 J 

Convergent : 

-1.041 

-,05848 

-.3639" 





G 

—ns 

-.4058 

-.02163 

-.1464 







176.6 

36.37 

55 . 60 

II 

- o 

36.37 

9.797 

12.06 

1 

55.60 

12.06 

21.81 


176.9 

36.39 

55.71 

1 

= 

36.39 

9.798 

12.06 


55.71 

12.06 

21.85 


197.4 

37.56 

62.83 

1 

^2 " 

37.56 

9.868 

12.46 


62.83 

12.46 

24.35 


166.4 

35.79 

52.08 


35.79 

9.762 

11.86 


52.08 

11.86 

20.58 

_ 



177.7 

36.43 


36.43 

9.801 


55.98 

12.08 


A 


c 


55.98 12.08 21.94 
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T 

Expected cost = x ^ jc 
Stability: 


Configuration 


0 ) 

1 ( 5 - 1 ) 

2 (B 2 ) 


3 (B 3 ) 


Stable 

-yes- 

yes 

yes 

no 


Interpretation : The system 

= [A + B . G ] 

— 1+1 — — 1— ns — t 

is stable for i = 0,1,2. 

Configuration state 2 is stabilized because the probability of the 
configuration state being 2 (^3 ^ larger than in Case i) . 



Example 6.2 Case iii) 


= .01 p = -98 

f r 

0 0 


p = .10 p = .90 

1 1 


7T = 


.8909 

.009172 

.09891 

,001010 


Convergent: 


G 

—ns 


-1.729 -.09453 -.6062 

-.3400 -.01858 -.1195 




210.6 

41.04 

67.28 


1 

^0 

= 

41.04 

10.76 

13.61 




67-28 

13.61 

26.29 




213.2 

41.14 

68.26 



= 

41.14 

10.75 

13.66 




68.26 

13.66 

26.66 






— * 




212.3 

41.09 

67.92 


1 


41.09 

10.75 

13.64 




67.92 

13.64 

26.53 




196.0 

40.19 

62.11 


^3 

= 

40.19 

10.70 

13.32 




62.11 

13.32 

24.47 





210.7 

40.99 

67.28 

So 

IT. S 
1 — 

r 

i 

40-99 

10.75 

13.60 




67.28 

13.60 

26.28 


& c 
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T 

Expected cost = x C x 
Stability: 


Configuration 

Stable 

o 

o 

-yes 

1 (Bj) 

yes 

2 (B^) 

yes 

3 (B ) 

Interpretation; The system 

X 4-j.T = [A + B . G ] X ^ 

— 1+1 — — 1 — ns — t 

is stable for i = 0,1,2. 

no 


Configuration state 2 is stabilized because the probability of the 
configuration state being 2 (B ^2 ^ larger than in Case i) - 
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6.4 Summary . 

In this Chapter, two applications of the non-switching gain method- 
ology to computer-aided design (CAD) were presented. The purpose of 
these examples was to demonstrate the usefulness of the non-switching 
gain methodology in the design process. CAD has two uses: First, it is 

used by the system designer in the evaluation and design of a system. 
Second, it is quite useful to the theorist. In this research, for 
example, without CAD techniques, a thorough knowledge of the methodologies 
presented in this report could not have been gained. The equations 
describing the switching and non-switching gain methodologies can be 
derived, but their meaning in a specific context cannot be determined 
theoretically. The purpose of this research was to study the inter- 
actions between system reliability and optimal control. The method- 
ologies presented in this report allow this study to proceed. The two 
Examples of this Chapter study two specific areas of interaction 
between system reliability and control. The door has now been opened to 
the answers to questions on reliable control system designs. Computer- 
aided design can provide the signposts to these- answers. 
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CHAPTER 7 
CRITIQUE 


7.1 Introduction . 

In this Chapter, the major results of the report will be summarized. 
In Chapters 3 and 4, the switching gain solution was developed and 
extended suboptimally to stochastic systems. In Chapter 5, the non- 
switching gain solution was developed. The problems associated with 
system stability, including definitions of what constitutes a stable 
system, and with the steady— state solutions to Problems A (Sections 3 
through 5) and B (Section 6) were studied in detail in Section 7. The 
equivalence of the two approaches to the non- switching gain solution is 
proved in Section 8. The existence of a robust steady-state linear 
feedback control system was studied in Section 9. 

In the following sections, each major result will be discussed; in 
Section 5, some suggestions for future directions in research will 
be made. 

7.2 The Switching Gain Solution . 

The switching gain solution was derived in Chapter 3 as a control 
methodology for linear system with quadratic cost criteria and variable 
actuator configurations. The resulting control law was to account for 
the failure, repair and reconfiguration of the actuators by switching 
the control gain on detection of a change in configuration. This type 
of control law is, from Chapter 1, Section 4, a class II reliable control 
methodology; an active (switching) controller is used with a passive 
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configuration design. 

7.2.1 Deterministic Optimal Solution - 

The switching gain solution of Chapter 2 is derived as the optimal 
solution for the discrete-time deterministic optimal control problem. 

It is the optimal control simply because the structure of the discrete- 
time system allows perfect observations of the system structure with 
one-step delay. Therefore, there is no need for the control law to 
have a dual effect; in fact, there can be no dual effect, since the 
control law does not affect the observation process, for almost all 
values of the control. 

A minor drawback to the switching gain solution is the computa- 
tional burden of iterating the Riccati-like equations (3.3.6), and solving 
for the optimal control using equation (3.3.7), backward in time for 
each time instant of the control interval, or until the steady-state 
solution is achieved, when one exists. Fortunately this computation is 
done off-line, and the various optimal gains are then stored for on-line 
use. On-line, the controller simply determines which structure the 
system was in at the previous time instant and chooses the corresponding 
(stored) gain- The control law is then a linear feedback control using ' 
that particular gain. 

7-2.2 Non-Extendability to Stochastic Systems - 

Unfortunately, the switching-gain solution does not extend optimally 
to systems where noise is present. When noise is present, it is no 
longer possible (in general) to determine exactly the previous value 
of the system structure. It was shown in Section 3 of Chapter 2 that 



176 


in such a case, the optimal control law exhibits a dual effect; i.e., 
the control law influences the measurement of the system structure. In 
a real-life situation, it is unlikely that a system with no internal 
noise will be found. Unfortunately, the optimal (dual) control law is, 
in practice, unsolvable due to the immense computer resources which are 
required. 

7.2.3 Suboptimal Extensions. 

Because of the dual control effect, the deterministic optimal 
solution is the only closed-form solution available. Thus, it is in 
our interest to look for suboptimal methodologies which extend the 
switching gain solution to the stochastic case. In Chapter 4, two of 
these methodologies were studied: Hypothesis testing and dual identi- 

fication. While hypothesis testing is a measurement strategy, dual 
identification modifies the control in order to guarantee a perfect 
observation of the system structure with the next measurement. Both 
methodologies are presented in their simplest form, since the problems 
of stochastic control of systems with variable structure are not within 
the scope of this research. Two comments are in order, however: 

First, at least in the form presented in Chapter 4, a dual identifica- 
tion algorithm is computationally intensive. Since it is an on-line 
algorithm, a significant computational capacity may be required in its 
implementation. Second, it is observed that the optimal stochastic 
control law, if it could be calculated, would rely on both estimation 
and dual control, the two concepts which are represented in Chapter 4 by 
hypothesis testing and dual identification, respectively. 
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In a suboptimal implementation using dual identification, the 
algorithm would most likely be used only at intervals; the implementa- 
tion would rely on an estimation algorithm for the remainder of the 
time. This scheme would attempt to minimize the degrading effect of 
dual identification on the state trajectory by using it only to guarantee 
that the estimation algorithm was tracking the system configuration 
properly. Thus, the system response would be roughly periodic, with 
the state being driven away from the origin in order to obtain an 
accurate estimate of the configuration, and decaying back toward zero 
between uses of the dual identification algorithm. 

This type of control strategy deserves some attention in future 
research activities. It is similar to the class of self-testing 
systems which perform diagnostic testing of their configurations 
at intervals. It is also, at present, the only methodology which takes 
advantage of the dual property of the control law in systems with 
variable, imperfectly observed, structure. 
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7.3 The Non-Switching Gain Solution . 

The non-switching gain solution of Chapter 5 was derived as an 
alternative to the switching gain solution of Chapter 3. Although 
the non-switching solution is, in general, suboptimal, the on-line 
complexity of the solution is less demanding than that of the switching 
gain solution. On-line, the non-switching gain solution has the same 
complexity as does the standard linear quadratic solution. Off-line, 
the computational requirements are equivalent to those of the sv/itching 
gain solution. 

7.3.1 The Necessary Conditions — Unsolvability. 

When the non-switching control problem is formulated as an 
equivalent deteanninistic control problem (Chapter 5, Section 4), the 
necessary conditions from the matrix minimum principle [Athans,41] 
yield a two-point boundary value problem which is not explicitly 
solvable; at the present time, the solution to this problem appears 
intractable. The necessary conditions are used, however, in conjunction 
with an equivalent problem (Chapter 5, Section 6), to prove some strong 
properties of the solution to the equivalent problem. 

7.3.2 The Equivalent Problem. 

The equivalent problem formulated in Section 6 of Chapter 5 has 
the advantage over the original formulation that a closed-form expression 
for the solution can be readily obtained. From the necessary conditions 
of Section 5 in Chapter 5 for the original formulation, it is shown that 
the accxomulated costs over the control interval for a specified gain 
sequence are identical for the two formulations. Prom this, in Section 8 
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of Chapter 5, it is shown that if the steady-state solutions to both 
problems exist, then they are identical- This is a major result, since 
the steady-state solution to the second formulation is calculable, 
while the solution to the first formulation is not. 

7.3.3 Existence of a Stabilizing Gain. 

Only one major result remains; one would hope that the steady-state 
solution to the second formulation exists if and only, is the steady-state 
solution to the first formulation exists. In Section 7 of Chapter 5, the 
meaning of "steady- state" is precisely defined for both problems- In 
order for the concept of a steady-state solution to be well-defined, an 
exact definition of stability must be given- Two definitions are present- 
ed. Stability is defined as the usual concept of mean-square stability. 

A definition of cost-stability is presented as the condition when the 
expected cost for the infinite horizon problem (unno 2 unalized by time) 
is bounded. It is proved that the solutions to the two formulations 
are equivalent in that one solution is cost- stabilizing if and only if 
the other is also. Cost stability is shown to imply mean-square 
stability; the reverse is not necessarily true. 

7.3.4 Problems with Convergence . 

There are two criticisms of the results of Chapter 5. First, 
although cost-stability is not implied by mean-square stability, it is 
possible that, for the specific form of the non-switching gain solution, 
the two definitions are equivalent. This is a minor point, in that the 
equivalence result is already very strong; it yields a procedure for 
the calculation of the steady-state solution to the two point boundary 
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value problem which converges if and only if that solution exists. 

Second, there is still a minor problem concerning the convergence 
of the non-switching gain solution. The equivalence theorems of 
Chapter 5 only require the solution to have a steady-state, which may 
be a limit cycle. A limit cycle is still copacetic, but it is harder 
to implement than one gain would be . Therefore , it is desired that 
conditions be found for which the possibility of a limit cycle is 
ruled out. 

Thus, two possible topics for future research are the examination 
of the exact relationship between cost-stability and mean-square stability 
for the non- switching solution and the determination of conditions for 
which the possibility of limit cycles as solutions is eliminated. 

7.3.5 Existence of a Robust Gain. 

A spin-off of the non-switching gain solution of Chapter 5 is 
the development of an algorithm which determines when a robust gain 
for a set of linear systems exists (Section 9) . A robust gain is a 
gain which stabilizes each mode of the system configuration regard- 
less of the configuration dynamics. This algorithm is developed by 
noting that the robustness problem can be reformulated as a non-switch- 
ing gain problem. Since the non- switching gain is, in the steady-state 
case, the solution to the first formulation (Section 4, Chapter 5) , and 
since it is stabilizing if and only if a stabilizing gain exists, then 
by the special structure of the robust formulation (Section 9) , the 
steady-state non-switching gain is robust when it exists, in addition, 
if the non-switching solution is not cost-stabilizing, then no robust 
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gain exists. This is a very important result; it is unfortunate that 
determination of existence of the robust gain reguires the solution 
of the non-switching gain problem. At present, however/ ho test bn -a 
system exists which determines when the non- switching gain solution 
is cost-stabilizing'. It is hoped that such a test will be developed in 
the future . 
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7.4 Computer-Aided Design . 

Chapter 6 demonstrates the usefulness of the non-switching 
gain calculations in computer-aided design (CAD) . These calculations 
provide the backbone for comparison studies on the relative system 
effectiveness of various designs. In the first example, it is demon- 
strated that the non-switching control methodology yields a numerical 
value based on the expected performance of a design configuration 
over the effect of the structural dynamics. This example demonstrates 
that relatively subtle qualities of an actuator can be used to rank 
various actuator configurations; in this case, the ranking depends 
on the manner in which the control affected the system state and is 
not obvious on a casual inspection of the configuration. 

The second example demonstrates the ability of the non-switching 
gain methodology to observe the trade-off between high reliability and 
high effectiveness in an actuator. Both qualities are desirable, but 
in this example, one actuator is highly reliable,, while the second 
actuator is not as reliable, but is highly effective in-that it incurs 
a much smaller cost in applying the same control effect to the system. 

The non-switching gain problem is solved for a range of actuator reli- 
abilities for the highly effective sensor, it is demonstrated that 
the trend exists to depend more heavily on the high reliability sensor 
as the reliability of the highly effective sensor decreases, even 
though the operation of the highly reliable sensor incurs more cost. 

Chapter 6 only touches upon the field of computer-aided design. 

There is much work to be done in this field, and the purpose of Chapter 6 
is only to establish the usefulness of the non-switching gain methodology 
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in the design process. In the future, the applicability of the non- 
switching gain methodology to CAD should be studied in great detail; 
in particular,, a comprehensive methodology for the application of the- 
techniques of Chapter 5 to CAD should be developed. This methodology 
should include a strong argument for the validity of using the non-switch- 
ing methodology in CAD. Specifically, research needs to be carried out 
on the relationship of the costs incurred by various design configurations; 
this is similar to justifying the use of the quadratic cost criterion 
in the linear quadratic regulator. In order to compare two designs, a 
valid basis of comparison, or cost index, must exist. The non-switching 
methodology is proposed as being a valid cost index for the class of 
systems for which it is applicable; this conjecture should be verified. 

In addition to the usefulness of the non-switching methodology, it 
has been mentioned previously that a valid definition for a reliable 
design is that the design is cost-s tabilizable . Since, for the deter- 
ministic control problem presented in Chapter 3, the switching gain 
solution is the optimal solution, the existence of the steady-state 
switching, gain solution is equivalent to the stabilizability of that 
design. Hence, the existence of the steady-state switching gain solution 
is necessary and sufficient to classify a design reliable . 

In theory, the computation of the steady-state switching gain 
solution can be used as a method in CAD for determining if a proposed 
design meets the minimum requirement of stabilizability. In practice., 
however., the proposed design will operate in a stochastic environment; 
therefore, the switching gain solution is not an absolute measure of the 
stabilizability of the design. In the future, research should be 
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concentrated on the development of the concept of stabilizability to 
more general stochastic systems than has been done previously. An 
example of work in this direction has been given with the Uncertainty 
Threshold Principle [Athans, et . al . , 37] , which is basically the deter- 
mination of conditions of stabilizability for a specific system with a 
specific type of control law. The work on the existence of the non- 
switching gain solution for a simple system (Chapter 2,- Section 7) 
is another example. It has been demonstrated in this research that the 
concepts of systems reliability and stabilizability are crucially 
interconnected. It is left to future research to determine more general 
conditions of reliability and stabilizability and to implement these 
conditions in computer algorithms which can be used by the designer. 
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7 . 5 Suggestions for Future Research . 

Several suggestions for future research have been presented in 
Sections 2,3 and 4 of this Chapter. In -this Section, a summary of these 
suggestions will be given. 

In Chapter 1, three classes of reliable control methodologies 
were given. These are 

I) Passive (Robust) Controller Design 

II) Active (Switching) Controller, Passive Configuration 
Design 

III) Active Controller, Active Configuration Design 

Of the methodologies presented in this report, the non-switching 
gain design is a class I methodology, and the switching gain design is 
a class II methodology. Class III methodologies are not represented 
in this report. This class is currently largely in the realm of 
"blue sky" theory. Unfortunately, there is as yet no adequate model 
of configuration dynamics which exhibits a state and control structure. 
Over the next ten years, one should see much research activity in the 
area of class III methodologies and their control structures. 

In class II methodologies, much effort should be concentrated on 
extensions, either optimal or suboptimal, of the switching class of 
control laws to stochastic systems. At present, most work, has been done 
in estimation theory, since the difficulties associated with dual 
control are widely recognized. The ability of a control law to perform 
diagnostic testing for changes in configuration has yet to be exploited 
theoretically, although many heuristic algorithms have been used, both 
in control systems and in the more established field of fault detection 
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and identification in digital systems. Dual control is a form of self- 
testing, and can be utilized as such, even if an optimal control is 
not known. The dual identification methodology of Chapter 4 is an 
example. This field requires a large effort, and should be rich in 
research opportunities. 

The class I methodologies are represented in this research by the 
non-switching gain solution. The work done in Chapter 5 on mean-square 
stability and cost-stability of solutions is not unique to this class of 
problems. Much remains to be done in the classification of what consti- 
tutes a stabilizable system, whether with respect to a non-switching 
control law or something more general. 

Since reliability can be defined as stabilizability with respect 
to some class of control laws, research into the stabilizability of 
dynamic configuration systems is the key issue in reliable control 
system designs. Much work, including this research, has been done on 
the assiomption that the system is stabilizable; however, little progress 
has been made in determining why a given design is stabilizable. 

Although iterative tests were developed in this report for determining 
stabilizability, a thorough understanding of the reason these tests 
either converge or fail to converge is lacking. Much work still must be 
done. With this should come a resolution of the problems with limit 
cycle steady-state solutions to the non-switching gain methodology. 

In Chapter 6, the usefulness of the non-switching gain solution in 
computer-aided design was demonstrated, CAD is a field unto itself; many 
opportunities exist for research in this area. Unfortunately, most 
research is application— specific . CAD is useful not only to the designer. 
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but also to the researcher. It is a powerful tool in the building of 
the concepts of reliable control systems design, and it should be 
developed in parallel with any future .research. 

7.6 Summary . 

In summary, the main purpose of this research was to establish a 
foundation in reliable control system design methodology which would 
provide the basic concept of a reliable control system. In achieving 
this goal, the linear quadratic variable actuator control problem was 
studied in some detail. Optimization problems were formulated which 
represented both system performance (in the quadratic performance index) 
and .system reliability (in the expectation of the performance index over 
all possible structural trajectories) . The optimal control law was 
solved analytically for the deterministic system; this was the switching 
gain solution. It was clearly illustrated by example in Chapter 2 that 
the switching gain control law could not be extended analytically to 
the control of stochastic, systems . This example demonstrated the dual 
effect of the control law; in general, the control law will influence 
the measurement accuracy optimally (in the sense of minimizing expected 
cost) when the control can influence the accuracy. 

Stochastic extensions to the switching gain methodology were proposed 
in Chapter 4. In particular, the dual identification algorithm is an 
illustration of the self-testing capacity of dual control laws. The 
study of the uses of the dual control effect in the design of reliable 
control systems is a promising research area of the future. 
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In Chapter 5, the non-switching gain solution was developed. This 
solution led to an algorithm for the determination of robust linear 
constant gain control laws for a set of linear systems with different 
actuator configurations. In addition, the resulting gains are optimal 
with respect to a given quadratic performance index and exist if and 
only if any robust gain exists. 

In conclusion, the unifying concept of this report is: What 

constitutes a reliable control system, or a reliable design? A major 
connection was established in this research between the concepts of 
reliability and stabilizability. Iterative procedures were developed 
for the determination of whether or not a given linear system of the 
type considered in this report is reliable, with respect to both class 
I and class II controllers; i.e., non- switching and switching gain 
controllers, respectively. 
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DEFINITIONS FROM MIL-STD-721B 
25 August 1966 


RELIABILITY 

The probability that an item will perform its intended function 
for a specified interval under stated conditions. 

AVAILABILITY 

A measure of the degree to which an item is in the operable and 
committable state at the start of the mission, when the mission is 
called for at an unknown (random) point in time. 

DEPENDABILITY 

A measure of the item operating condition at one or more points 
during the mission, including the effects of Reliability, Maintain 
ability and Survivability, given the item condition (s) at the start 
of the mission. It may be stated as the probability that an item will 
(a) enter or occupy any one of its required operational modes during a 
specific mission, (b) perform the functions associated with those 
operational modes. 

CAPABILITY 

A measure of the ability of an item to achieve mission objec- 
tives given the conditions during the mission. 

OPERABLE 

The state of being able to perform the intended function. 
MAINTAINABILITY 

A characteristic of design and installation which is expressed 
as the probability that an item will be retained in or restored to a 
specific condition within a given period of time, when the main- 
tenance is performed in accordance with prescribed procedures and 
resources. 

SURVIVABILITY 

The measure of the degree to which an item will withstand hostile 
man-made environment and not suffer abortive impairment of its 
ability to accomplish its designated mission. 
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A2.1 Exact Optimal Solution for Deterministic Case , Chapter 2, 

Section 2 . 

From (2.2.7) and using dynamic programming, we wish to minimize 

- 2 2 

V{x , k(t-l), u , t) = E(qx, + ru 
t • t t t 

+ V (ax +b, , . u ,k(t),t+l[x ) (A2.1.1) 

t iZ. vt } t t 


where V (*,k(t), t+1) represents the minimum cost-to-go, given 
k(t) at time t+1- 

This minimization can be carried out because is known exactly 
at time t, and therefore ^is knownexactly by equation (2.2,10). 

The control u^ is computed from 


0 = 


* 1 \ 

ir^ V (a x^+, k=l,t+l) j 


2 * 

+ ru^+TT^ V (a x^+ b u^, k=0,t+l) 
t 


(A2.1.2) 


and the assumption that 


(A2.1.3) 


resulting in equation (2.2.8). Equations (2.2.12) and (2.2.13) are 
then obtained by substitution of (2,2.8) into (A2.1.1) j these 
equations validate assumption (A2.1.3) by induction. 


5 
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a 2 . 2 Exact Optimal Solution for Stochastic Case , T=Q , 1 , 2 = 

(1-d example) . 

The formulation is the same as in A2.1, except the system is 
now represented by 

\(t) “t 

is white noise with zero mean, variance 5, and probability dis- 
tribution p(5), which is uncorrelated with any other variable. To 


illustrate the complexity of the solution, the time set is chosen as 
{ 0 , 1 , 2 } . The problem is to find u^ and such that 

V(x ,0) = E(J) = eI ^ (x^qt u^r) tx^q |x ,TL (A2.2. 

L t =0 


r) tx^q 0 


is minimized. Let V denote the minimum value of V. Assume 

u_j. = ^^t^ 

where 4> ^ is a mapping from the information at time t (2^) into the 
control space. 


\ ■'^”0' *0' “0 “t-i' 


then 

A (2*2^ 1 

V {Xq, 0) = min E x^^q + u^r + V (x^,l)| 

by dynamic programming. Also 


V (x^,l) = min E x^ q + u^ r + V (X2^2)j | 

Ui=<jJi(2i) ^ ^ 

* 2 

But V (A2.2.6) becomes 

* (2 2 2 ,) 

V (x^,l) = min E x^ q+ u^ r + x^ q| Z J 

Ui=<()i(Zi) ^ > 




jxjqt ujr+ (3*1+ Uj +5 q|Zj^j 


(A2.2.8) 
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now 


' h = -^2:0 

(A2.2.8) = min / 


+ u^r 


+ E 


^ iT^(l|l) (ax^+ b^u^+ 5^) 
1=0 


(A2.2.9) 


where 7T^(ljl) is the probability that = i, given Bringing the 

expectation inside the sum. 


(A2.2.9) = min 


(Zj_l 


2 2 
x^q + u^r 


1^ 


(ill) (a^x^ + b?u^ + c + 2ab.X-U-)q, 
‘ 1 1 1 1 1 1 1 


(A2.2.10) 


Differentiating (A2.2.10) w.r.t. and setting the result equal to 


zero: 


0 = 2ru, + 7 IT , (1 ] 1) (2b .u, + 2ab.x.)q 

X X XX X X 

i=0 


(A2,2.11) 


or 




1 I 

^ IT dll) b qa 


r + 


iT.(i|i)b^ 


X, 


(A2.2.12) 


Substituting (A2.2.12) back into (A2.2.10), define and as 

= 5q . (A2.2.13) 


= (a + l)q 


- a 



(A2.2.14) 
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V 




(A2.2.15) 


A few remarks must be made about the probability distribution over k , 

t 

given or 


Notation: 

■n'^{tjt) = probability that = i, given the available information Z^. 

TT^(tIt+l) = probability that k^ = i, given the available information 
^t+1’ 


Prom the Markov property, 
7T(tit) = P7T (t-l|t) 


(A2.2.16) 


Equation (A2.2.16) is the propagation equation for the distribution TT^. 
The foinn of the update equation is given and proved in the following 
lemma: 


Lemma A2.1: 


7r.(t|t+l) = 




(A2.2.17) 


Proof : 


Note that 

p(xt^l-axt-b.ut) = p(x^^^lz^,u^,k(t}=i) 

where u^ is not a random variable. Also, 
7T^(t|t) = p(k(t)=i[z^) 

p{k(t)— i,TTQ , Xq r • * • r 
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then ( A2 .2.17) becomes : 


P<k(t)=i|z^^^) = 


p(x^_l^^|z^,u^,k(t)=i)p(k{t)=i|z^) 


»(== t+il V 


which is Bayes rule. 


Q.E.r 


Returning to equation (A2.2.5), and substituting (A2.2.15)/ 

E + ujr + xjs + T^|Zq| (A2.: 


V (x , 0) = min 

“ o " 1 “ o « o > 


= min 


u =d)„ (Z„) 

0 ^0 0 


2 2 


+ X. 


q(l+a^) 




' 2 2 
q a 


r + 


TT. (l|l)b1^ 
1 ' 1 


0i 


(A2.: 


= min 


u (Z- ) 
0 ^0 0 


X^q + ugr + Hq 



i-TidlDb^: 
1=0 J 


■ 2 2 
q a 


r + 


■ k-k- 
1 0 


where 


i,ii(i|i)bj 

Li = 0 x 


IT. 


kQ ,0 


(A2.: 


T^dlD = i 


3^ 


p (x, -ax„-b .u„)lT . ^ 

1 0 3 0 ] ,0 


p(x^-ax^-b.UQ)TT.^Q 


(A2,: 


.18). 


.19) 


. 20 ) 


. 21 ) 
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Equation (A2.2.21) is a combination of equations (A2.2-16) and (A2.2.17). 
Equation {A2.2.20) can only be solved numerically (in general) ; this 
requires a numerical minimization of a f.unction .the .computation of which 
requires four numerical integrations — a difficult task. 
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A2 . 3 Exact Solution of Stochastic Case Over T = 0, 1, 2 = T j. 
for a Specific Form of p (£) , Chapter ^ S ectioh 2.3.1 . 

Assume , for the problem in A2 . 2 , that 

2 ^^ , for -/3c £ ? £ 


P(5) = 


0 , otherwise 


Suppose [uqI > 0 is large enough such that 

p(.(bj^ 5 q) = 0, i / and £ [-i^,/3c] 


Then 


(A2.3-1) 


VirdlDb. = 




• . . P (x^ -ax -b .u^) ir . ^ 

^=0 ^ ° ^ ° ^ \ b^ (A2.3.2) 

P‘=‘r“o’bk“o?\,0 



0 


(A2.3.3) 


(A2.3.4) 


Similarly , 


TT. {ijDb^ = P-1, 

1 ‘ 1 ik„ 1 


0 


(A2.3.5) 


Then, from equation (A2.2.14), 


Si(ko) = (a^+l)q - a^ 


1=0 0 , 


r + 


% bj 


(A2.3.6) 
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2 

Uq^h- r;q 


+ b u + E ) 



dp(Eo), 


■^0 


TT, 


fco-o 


(A2.3.7) 



Differentiating with respect to u^, and noting that does not depend 


on Uq, 




X ^ X ^ 

~ ^ j ^It 0 ^ ^ If ^^^]f ^0 (A2.3. 

° k^O k“^0 Vo ^0 ° ^0 ° ^ 


9) 


Then 
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* 



1 ± 

£ \ ,0 £ p 

^0 k^O 


, , b- S, a 

Vi ‘'o ^ 


x_ 


— 2 

1^0 '" O '” 1^0 10 ^0 ^ 


■^\E„v.o X„ 


(A2.3.10) 


This solution is valid only when [ j > 0 is large enough such that 
p((bj^ -b^)u^ + ?q)- 0, arid £ [- /3^, /35] . Thus, 

1 -b^) Uq+ CqI >v^, [-v^ , /^l (A2.3.11) 


must be satisfied. 

i) Assume {b, -b.) >0. Then (A2.3.11) is satisfied if 

0 O 

(b, -b.) - /3f > /iE (A2.3.12) 

^0 ^ “o 

or 

(b, -b.) > (A2.3.13) 

'‘o ^ “o 

ii) Assume (b, -b.) <0, Then (A2.3.11) is satisfied if 

X u„ 

0 0 

(b, -b.) + /3H < ~/3" (A2.3.14) 

''o ^ “o 

or 


(b, -b.) < ~2y^ (A2.3.1S) 

"o ^ "o 

■* 

Therefore, u^ must satisfy 

I (b, -b.)u*| > 2t/^ . (A2.3.16) 

' ko X o' 

for (A2.3.10) to hold. 

'k 

Notice also that when (A2.3.10) is the optimal solution, u^ is 
identical to the deterministic solution. 
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A2-4 Existence of Steady-State Solution for 1-d Example . 

From Chapter 2, Section 2.2, the coupled Riccati equations for 
Sq and are 




^21 r ■ 




S = „ t ^ P22'°-^>^l.t4.l’ 




+ p /, ^g22<»'^>^.fH 

\ b[P * H. P2,(l/b')S^_^^,J 


\ = i: 


(A2.4.3) 
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S 


Q,t 


(A2.4.4) 


0,t+l 


Dividing both sides of equations (A2.4.1) and (A2.4.2) by t+l' 
manipulating terms, and using equations (A2.4.3) and (A2-4.4) yields: 

,2 


^t S 


r[Pj^j_ab +P2i<a/b)hj^^^l 


0,t+l 


’o,t+i P21 


+ P^l I a - 


b[p^j^ab+ 




+ £>21 






t+1 


(A2.4.5) 




_ q 


o,t+i 


+ 


1 bEPj^ab + P22<»/b>bt+ll 


■0,t+l + Pjjb^ + Pjjd/b^)!!^^,)^ 


+ *■ 12 '“ 


b[Pi2ab + p^2 (a/b)h^_l_j^] 

^^'^^0,t+l^ ^^12^ ■^^22^^'^^ ^\+l 


t+1" 
2 


+ P22«^ 


^12^^ ^22^^'^^\+l 


E’ 12 '’^ ■" ® 22 **t+l^^’ 


t+1 


(A2.4.6) 
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Assume S„ , S, -><» as t and h ->■ h, T F . 

0,t l,t t t: 


/ b[p^^ab + P2^{a/b)h] \ 2 

r = p I a - 2 2 1 

V ^ 11 ^' •*• / 


/ P22^(a/b)h \ 2 

■*■ ^21 I ^ ~ 2 2 1 ^ 

b[p^^b +P2ih/b]/ 


/ b[p^2^^ \2 

~ ^12 ( ^ ~ 2 2 I 

\ Pl 2 *» / 


/ Pl2^ P22(a/h)h ^2 

* ®22 r ■ ~' ";2 TTJ , I ** 


P = 


^11 Pl 2 


^21 ^22 


l] _ pi ^-I >2 

J L ^‘^1 ^ 2 . 


Then 


r= p la - 


b[p^ab + (1-p^) (a/b)h] \2 
p^b^ + (a/b^)h / 


/ P,ab + (1-p ) (aA>)h \2 

+ (1-P, ) 1 a - 2 2 I ^ 

\ b[p^b^ + (l-p^)h/b ]/ 


/ b[(l-p )ab + P 2 (a/b)h] \2 

hr = {1-p ) ( a - 5 5 

\ (l-p2)b^ + p^d/b )h3 / 


/ (1-p )ab 4- p (a/b)h 

•*• P2 “ 2 2 

\ b[(l-p2)b^ + p^d/b )h] 


(A2.4.7) 


(A2.4.8) 


(A2.4.9) 


(A2.4.10) 


(A2.4.11} 
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Solving for h and F from equations (A2.4.10) and (a2. 4.11) , if F > 1, 
then there exists no steady-state solution. 
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A3 - 1 Proof of Theorem 

Assume for k^S,. Then (B - B^)u^_^ = 0, 

which implies u is in the null space of B - B» ,.N(B - ^q). 

Now, dimension ^j^)) ^ because the distinct. 

Therefore, 


dimension {U N(B -B„))<m 
k,£ 


m 


Therefore the set U N(B -B. ) has measure zero in R . 

k,£ ^ 


(A3. 1,1) 
Q.E.D. 



A3.2 Optimal Solution for Deterministic Problem. 


For the system 


-t+1 -k(t) -t 


-k(t) ^ ^ -k ^k=0 




e R 


.L+l 


(A3.2..1). 
(A3. 2. 2) 
(A3. 2.3) 


where tt . = probability of B . at time t. 

X / 1 i 

Assume that 


1) JC ^ is observed exactly 

2) then B, ,, chanqes to B. . 

k(t-l) ^ — k(t) 

3) then u ^ is applied 

From dynamic programming, the optimal cost-to-go at time t is given 
by 

E. 


V (jt^ ,k(t-l),t) = min 


T T 

— t^— t — t t 




+ V ,k(t),t+l) x^ 


(A3. 2. 4) 


Assume 


V*(x ,k(t-l),t) = x^S, ^x_^ 

— t — t— k,t — t 


(A3. 2. 5) 


Then 




min 




T T 



IUqII 

< foj. all 

1 -11 F I 
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(A5.1.10) 

Q.E.D 
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A5.1 Proof of Theorem 1, Chapter 5. 


w 


T-1 




tr[^^ {Q+^ 


and J^<B. Since Q + ^ ^ constant for all 

iim tr[^ 1 = 0 

t-jco 


which is exactly Definition 1. 

(=^) 

From equation (5-4.6), note that 


— 1 , t+1 i=0 — 1 / 1 1=0 


where F(-) is linear in (E . ^ }. 

— 1 , t 1=0 


Since 


lim tr[^ ] - 0 

t-J-oo 

for any choice of , !! ^ !i bounded and ||p |[ < 1. 

3- Zq 

Then 

T-1 

- J = V - tr[E^ (Q + R G ) ] + - tr[E^Q] 
n T n — t — n 


< V ^ tr[E 1 ^ tr[Q + g'^RG] + ^ 

— n — t n — n 

t=o 


tr[2^] ^r[Q] 


W^f 11 sji liQ t gVII -^Ikf iisji IIqII 


I|f|| ^ ( !! E i| 11 Q + g'^RgII ) 


1 - 


1 - 


1 q II II Q' + 1 RgII 


(A5.1.1) 

t, this implies 
(A5.1.2) 

(A5-1.3) 

(A5-1.4) 
(Otherwise , 

(A5.1.5) 
(a5. 1.6) 
(A5.1.7) 

(AS. 1.8) 


T+1 


(A5-1.9) 
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and 


(A3. 2. 6) = min 


UIJ.11 -V 


( ^ 

/ \r 


x.Qx. + u.Ru. 
— t— — t — t 1 


y.p 


1=L 


ik 


T T T T 

x^AS. Ax^+u^B.S. 

— t i,t+l 1 — t ~i — 1, t+1 — 1 — t 


T T T T 

+ x^AS. + u^B.S. _ 

— t i,t+l— 1— t — t— X— i,t+l 


Ax 


t jj (A3. 2. 


7) 


Differentiating the r.h.s. of (A3. 2. 7) w.r.t. and setting equal 


to zero: 


r T T 

= 2Ru^ + >,p., 2B.S. + 2B.S. A 

1 xk —X — X,t+1”-X — t — X— X,t+1 — 


(A3. 2. 8) 


or 


u 


k (t-1) , t 


T 

^ ^4Pik®i"i,ttl®i 


-1 


i: 


p„B,S. .,,Ax. 
xk— X — X , t+1 1 


(A3. 2. 9) 


is the optimal , given k(t-l). 


Since no noise is present in the system, k(t-l) is obtained from 

X and X , , , along with u , , as 

— t — 1-1 ^ — 1-1 

k (t-1) = i iff x^ = Ax^^ +B.U. T (A3 . 2 . 10) 

— t 1-1 — X— 1-1 

Substituting (A3. 2. 9) into (A3. 2. 7), and eliminating x^ because the 
equation must be true for all x^ and the matrix equation is symmetric. 
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on simplification we obtain 



which verifies assumption (A3. 2. 5) by induction, along with the initial 
condition 

^ (A3. 2. 12) 
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A3. 3 Proof of Lemma 

Consider the optimization of the cost-to~go given k (t-1) at time 
t with final time T. This optimal cost-to-go is simply 


(x ^ ,k(t-l) ,t) 


(A3. 3.1) 


where T denotes the final time. For the process with final time T+1, 
the optimal cost-to~go is 


/k(t-l) ,t) 


= E 




(A3. 3. 2) 


Since this optimal sequence is not necessarily optimal for the problem 
with final time T, it must not incur less cost over’ {t, . . . , t). 

^T+1 ^-t 




I ^ ip 

^ Qx 


(A3. 3. 3) 


Since the expectation term of equation (A3. 3, 3) is non-negative. 




(A3. 3. 4) 


Now , note that 

V* (x^,k(t-l),t) = Xt^i^t it 

f T 


(A3. 3. 5) 


and that equation (3.3.6) depends only on the number of iterations 


(T-t) for the calculatxon of S . , and therefore, 

T 

"k "k 


(A3. 3. 6) 
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Therefore, {s ^ is an increasing sequence in that 


S . ^ - S . > 0 (A3. 3. 7) 

— x,t~l ~i,t — 

St 

Since, by hypothesis, V is bounded over t, the S . ^ converge. 

T — l,t: 

Q.E.D. 
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A. 5.2 Proof of Remark on Theorem 1, Chapter 5. 


w. 


and 


^ tr[E^£] 




Since Q > 0 
T 

t=0 

Therefore 


tr[^^ ] is bounded. 


tr[Z^ ] ->-0 as t->«» 
— t * 


The reverse implication is shown to be false by example* 


(A5-2.1-) 


(A5.2.2) 


(A5.2.3) 


(A5.2.4) 


Example 1 : Consider 


^t+1 “t 


"t 


Then 


r 2, 

E[x^] 


S ^ 
t+1 0 


0 


but 


^ 2 

gE[x^] 



E 

t=0 


t+1 


(A5.2.5) 

(A5.2.6) 

(A5.2.7) 


(AS. 2. 8) 


* 


Example 1 is provided by Dr. D. Castanon of ESL. 
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A5 . 3 Proof of Theorem 2 , Chapter 5 . 

Let I = {o, 1, 2, l} 

and 

£-”(1) = {(k(0), k(l)r k{i)'e i} 

CO 

Define the function y on the cylinder sets of £. (I) 

K = { (k(0) , k(l) , . . - ) i k(i) fixed for i < T } 
for arbitrary" T by 

y(k) = \(o)q ^k(l)k(0) * ^k(2)k{l)* “* * %(T)k(T-l) 
where ]T ^ is the initial probability distribution over I and P = j ) 
is the stochastic matrix of transition probabilities for the Markov 
chain. By a theorem of Andersen and Jessen [Loeve, p.91,42], this 
function defines a measure, y , on the O' -algebra of £^(1) generated by 

00 oo 

the cylinder sets, cr(j6. {!))- Since (I)) p 1, from the definition 

oo • * * ' ' 

of y on the cylinder sets of H (I)-, ' 

y: aa”(I)) -5- [0,1] (A5.3.5) 

is a probability measure, and since y extends uniquely from the cylinder 

OO 

sets, it is the probability of occurance of elements of '0 (£ (I)). 


-(A5.3.1) 


{A5.3.2) 


(AS. 3. 3) 


(AS. 3. 4) 


Let 

J^(x) : [0,«] (A5.3.-6) 

+ x^Qx,j, (AS. 3. 7) 

where 

X = A X + , u (AS. 3. 8) 

-t+1 1 -t 
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H-t "" ^t-t (A5.3.9) 

x= (k(0), kdl, k(2),...) (AS. 3. 10} 

and let 

J = lim J (A5.3.11) 

since is constant on the cylinder sets with fixed sequences of 
length T+1, is measurable. (There are a finite number of such 
sets.) By Theorem A of [Halmos, p.84,10],J is measurable with respect 
to y. 

CO 

J(*)(x) : 5, (1)^ [0,«>] (AS. 3. 12) 

Let 

= {x eX. (I) I J(x) (x) < <» for X £ r” } (AS. 3. 13) 

and 

X^ == ^'^(I) - Xj^ (AS. 3. 14) 

Then and X 2 measurable subsets of fi,°”(I) , and therefore 

E[J] < 00 ° (AS. 3, IS) 

because J (x) is a non-negative function on R^. 

But 

E^[E[J]] (AS. 3. 16) 

from equation (5.7.14), and by hypothesis, r.h.s (A5.3.16) is finite. 
Therefore, any trajectory x is an element of with probability 1, 
and has finite cost. 

Therefore, cost-stabilizes (S.3.1) with probability 1. Q.E.D. 
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A?. 4. Proof of Theorem 3, Chapter 

Notation ; In the proof, the sequences ns ^^0 

t 

■k 

be referred to by G and G respectively. 

— — ns 


Proof ; 

I) f=>) Suppose G is cost- stabilizing. Then J(G ) < <», 
' • ' — ns — ns 

* * * 

But G minimizes J. Therefore, J(^ )< J (^^^ ) => J }< °°. 

k 

Thus, G is cost— stabilizing- 


* * 

I-I) (<^) Suppose G is cost-stabilizing. Then J(^ )< where 


J(G ) = lim J (G ) 
T-x» 


CAS -4.1)' 


Since E [J (G) ] = J^(G), 
X ns^ - T - 


* * * 
J(G ) = lim E [J (G )] = E IJ (G ) ] 
~ X ns„ — X ns — 


(AS. 4. 2) 


which implies 


ns — 


(AS. 4. 3) 


Since G minimizes J , then 
—ns ns 

J (G ) < J (G*) < «> 
ns —ns — ns — 


(AS. 4. 4) 


and, since E [J ] = for all T, for fixed G, 
X ns„ T “ 

T 

j(G ) < ”. 

—ns 


(AS. 4. 5) 


which implies that G is stabilizing. 

— ns 


Q-E.D. 
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For the control interval starting at time 0 and ending at time T, 


the expected cost for the optimal control is 

T 


(A5 .5.1) 


from equation (5.5.8) , where the subscript T refers to the endpoint 
of the control interval. Similarly ,, for the same process ending at 
T+1, the optimal expected cost is 




(A5.5.2) 


= E 


T 


T * T X 

(Q + G^CIM-I) R G^( T+1) 


+ =iT+ia2.I+ll^0' ^0 


= E 


fp * m * 

x^ {£ + G^ (T+1) R G^(T+1) 


x"^QX l2„,Tr„ 
— T — — T * — 0 —0 

_ 


(A5.5.3) 


+ E [ (G* (Ttl)'^ R G * (T+1) X^ + ^T+lS. 21 t+iI ^0 '-0 ^ 

(A5.5.4) 

The first expectation of equation (A5.5.4) is the cost corresponding 

* 

to the interval [0,T] , and must be greater than or equal to the 


second term is positive. Therefore, 

* * 

J , > J 
T+1 — T 


(A5.5.5) 


* * 

Since is bounded by hypothesis for all T, there exists a J such 

that 


= J 


(A5.5.6) 


* 

lim 

T 


Q,E,D, 
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a 5 . 6 Proof of Lemma chapter 
By direct computation, 

(G) = J„{G) + Efx'^ g\ G X + 

T+1 — 0?— ' — T ^ — qj — T+1 

and since the expectation is positive. 


Q X 1 
- - T+1 ■' 




since J^(G) is hounded, it converges. 


(A5.'6.1) 


(A5.6.2) 

Q . E . P . 
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A5 . 7 Proof of Theorem chapter 

converges to the steady- state value 
, and therefore , 


A) G 
— ns 


-G , because G 
—opt —ns, 


t ““t 

which minimizes the infinite-time horizon cost J 


ns 


ss 


by the argument given above, also minimizes equation (5.8.9) . 


B) Given e>0, a T>0 can be chosen which guarantees j| G^- ^ 

is ic 

11^ . - E . |[ < e and jl IT - it || <e , for all t > T. 

X / 1 X u 

f * 1 00 

Then, by the Principle of Optimality, the sequence i G^ ^t~T 

minimizes the infinite-horizon cost-to-go at time T. Consider the 

problem min J (G) for initial condition E . , ir , which has a solution 
^ ss - -1 - 

G independent of E . . In the limit as £-+0, the sequence l G, L . 
—ns —1 t^t=T(e) 


approaches the constant sequence of gains G . Suppose } *^^0 ^ 


T(£) ' 


the optimal cost-to-go, satisfies 

* * 

V^. < J - 6 

T(e) — ss 


(A5.7.1) 


Then the sequence of constant gains G would yield a strictly lower 


cost J (G ) 
ss — 


★ * 

J (G ) < J 
ss — ss 


(A5.7.2) 


since V . approaches the optimal cost-to-go, given the constant 
T (£ ) 

* 

sequence of gains ^ , in the limit, which is the solution to the 

equivalent problem min J (G) for initial conditions ^ , 2L- 

G ■ ^ 

Therefore 


* 

G = G 
— —ns 


(A5.7.3) 


Q.E.D 
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AM FCKTBAN 
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SUBRDUriNE AM(NAA,m,m,NQ,NR,NG,NS,KRA,N,M,iCCN,A,B,R,Q,P, 

1 3ST, E, S, SB, U, V,W, X, Y, PR, FZ , QWRM, lAD, lADINV, BS8,WCRK, IPVT, lEND, 

2 IPRT) 


** ***PARAM ETERS : 

INTE(ER l^,m,]NB,NQ,M,N3,NS,My\,N,M,K:CN,IPVT(N) ,IEND, mT 
DOUBIE PRECBI31SI BSB (NS,NAA,ICCN) , X CMA, N) ,EAD O.’^RA, N) ,iADINV (SRA,N) 
DOUBLE PRECISION E (KCQ^) ,SBT(NS,N) ,A{NA,N?VA) ,B(t©,lSAA,K:CN) 

DOUBIE fREClSION Q (NQ,N) ,R(NR,M) ,P a-^A, ICavi) ,S (NS,im, ICCN) 

DOUBIE PRECISION SB (NS,l«yi, ICCN) ,U(NA,N) ,V (NA, N) ,W{NA,N) ,Y (NA, N) 
DOUBIE IRECKION PR(N) ,WCRK(N) ,EZ{N) ,QSIORM(NG,NAA,ICON) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*****LOCAL VARIABLES: 

DOUBLE PRECISION COJD 

INTECER KIN, four, I, K, KMl, KK, J, JEND, L, KP, K41, ICTMl, IMl 
lOTECER ICOUNT 

IB ROUTINES CALLED; 

, M AED, M UNE Q, OPMAIB , MM UL , MXA IE , MATIO , EIGVA L, WEICHT , TRNATA 


*****PURFOSE: 

IHIS DOUBLE PRECISION StBRQUriNE CCMPUTES THE STEADZ-STATE OPTIMAL 
SOLUTION AND THE CCRRES PONDING OPTIMAL GAINS PCR. THE PRCBIEM 
DESCRIBED IN THE PUBLICATION: ' ON THE RELATIONSHIP BETWEEN 
RELIABILITY AND LINEAR QUADRATIC OPTIMAL CONTROL ' 

BY J. DOUGLAS BIREWELL AND M. ATHANS. 

(EQUATIONS (29) AND (30)). 


*****PARAMETER DESCRIPTION; 
CN INPUT; 

PAA 


NA,PB,NQ,NR, 

NG,NS,ISRA 


N 

M 

R3CN 

A 


THE SECQnID DIMLNSION OF THE ARRAYS S,SB,GNORM, 
BSB,B AS EECIARED IN THE CALLING PROGRAM 
DIMENSION STATEMENT; 

THE FIRST DIMENSION OF THE ARRA^ 

A (AND P, X,U,V,W,Y),B(AND BSB) ,Q,R,GNORM, 

S (AND SB,33T) ,RAD (AND RADINV) RESPECTIVELY 
AS EECIARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT ; 

THE NUMBER OF STATES ; 

THE NUT4BER OF CBSERmTIONS ; 

THE NUMBER OF CCNP IGURATIONS ; 

N BY N SSTEH MATRIX; 


ORIGINAL PAGE IS 
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AIM FORTRAN 


C 

C 

C 

■G 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

V 

c 

c 

c 

c 

c 

c 

c 

c 

c 


B 

R 

Q 

p 

E 

OUTHJT: 

m, rz 

U,V,W,S8T, 

X,Y 

s 

SB,BSB 

GKORI'l 

RAD^RADINV 

wceK 

IPVT 

lEND 

IPRT 


N BY M BY KCCN SET OF INFOT MATRICES ; 

M BYM CCNTROL WEIGHTING MATRIX;- 

N BY N STATE WEIGHTING MATRIX; 

KCa'I BY KCCN PRCBABILITY MATRIX; 

\ECTDR OF lENGTH KCCN CONTAINING THE NORMALIZED 
EIGENVECTOR OF P CCRRESIONDING TO THE EIGENVALUE 
CNE; 


SCRATCH VECTORS OF LENGTH N ;, 
N BY N SCRATCH ARRAS ; 


N BY N BY KCCN SET OF SOLUTIONS; 

N BY N BY KCCN SCRATCH ARRAS ; 

N BY M BY KCCN ARRAY WHICH WHL CONTAIN THE 
(AIN MATRICES FCR THE NORI4AL UNEAR QUAIRATIC 
(ALBS IAN PROBLEM; 

N BY N SCRATCH ARRAS ; 

SCRATCH VECTOR OF lENGPH N; 

SCRATCH VECTOR OF LENGTH N; 

NUMBER OF ITERATIONS USED IN SOLVING BOTH THE 
LINEAR QUAERATIC (ALBS IAN PRCeiEM AND THE 
HK»IEM DESCRIBED ABOVE ; 

FIRST ITERATION AT WHICH THE SOLUTIONS- W ILL BE 
HUNTED; 


CCMMCN/TNOU/KIN, K5UT 
ICOUL^T = 0 
DO 215 KK=1,IC® 

DO 4 J=1,N 
DO 3 1=1, N 

3 Y(I,J)=0.0O0 

4 Y(J,J)= 1.030 
DO 210 K=1,IEND 

CAIL MCF(«A,NB,NV,N,M,Y,B(1,1,KK) ,U,WCRK) 
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CAIL MADD a'lA, U, R, U) ^ 

DO 14 J=1,M QUALIIY 

DO 13 1=1, M 

13 V(I,J)= 0.030 

14 V(J,J)= 1.2D0 

CALL MLIKEQa^A,m,M,M, U, V,CCND, IPVT,WCRK) 

CALL TRNATB(NA,1«,N,M,B(1, 1,KK) ,X) 

CALL MMUL (NA,NA, N,M, N, X, Y, U) 

CALL MMUL (NA,m, m, N,M, N, U, A, X) 

CALL MCF(NA,NA,m,M,N, V,X,RAD,WCRK) 

CALL MSCAIE (LIRA, K, N, -1. 030, RAD) 

CALL MC?(NA,m,m,N,N, Y,A, U,WCRK) 

CALL MAID (NA,M, m, N,N, U, Q, U) 

CALL MAED m, NBA, m, M, N, U, RAD, Y) 

210 CCNTINUE 

KKMl = KK - 1 
WRITE (KOUr, 44441) 

WRITE (KOUr, 44442) KKMl 
CALL MATIO (NA,N, N, Y, 3) 

CALL MMUL{NG,NA,N^,N,M,M,V,X,GNORM(l,l„KK)) 

CALL MSCAIE (NG,M,N,-1. 030,QIORM{1, 1,KK)) 

CALL MT4UL(NB,1«3,NA,N,N,M,B(1,1,KK),(3I0RM{1,1, 1) ,V) 

WRITE (ROUT, 6000 ) 

CALL MATIO (NG,M, N, QK3RM (1, 1, KK) , 3) 

CALL MAED0SIA,NA, KA,N,N,V,A,V) 

WRITE (KOOr, 44443) 

CALL MATIO (NA, N, N, V, 3) 

CALL EIGVAL(NA,N, V,V,IR,K,WQRK,IPVr) 

215 CONTINUE 
JEND= 1 

WRITE (KOITT, 8000 ) 

CALL MATIO (NA, iCCN, iCOT, P, 3) 

DO 35 K=1,ICCN 
DO 30 J=1,N 
DO 40 1=1, N 

S (I, J,K)= 0.030 
40 CONTINUE 
30 S (J, J,K)= 1. 030 
35 CCNTINUE 

START ITERATION TO CALCULATE S (1 ) ,S (2) ,. . . S (K) ,Q3PT 

CAECULATE SB 
1 CONTINUE 
DO 50 k=i,h:cn 

CAIL MMUL(NS,]®,ffi,M,N,N,S(l,l,K) ,B(1,1,K) ,SB{1,1,K)) 

50 CONTINUE 

CAIL WEIGHT 0J3S,NAA,ICaSi,NS,N,M, E,SB, ST) 

CALCULATE RADICAL 
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DO 55 k=i,k:cn 

CAIL M ff (NS, IB , N, M, S (1 , 1, K) ,B ( 1 r If K) fBSB ( 1 , 1, K) ,WC»K) 

55 CONTINUE 

CAIL WEIGHT. (NB,.tC^,-ICCN,NRA,M,M,.E,BS8,-IftD) 

CALL MAED ^RA,NR,m,M,M,mD,E, U) 

DO 54 J=1,M 
DO 53 1=1 ,M 

53 RADINV{I,J)= 0.ffi0 

54 RADINV(J,J)= 1.0D0 

CALL MLINEQ(NA,NRA,M,M,U,RAEa;NV,CCND,IPVT,WORK) 

C 

c CArCULATE NEW SI, 1 = 1 , 2, ICON 

100 DO 1000 k=i,k:cn 

CALL MMULa'lS,NRA,NA,M,N,M, SBT,FADINV,U) 

CALL WEIGHT (NS,NAA, ICON, m, N,M, P (1, K) ,SB,V) 

CALL TRNA'IB{NA,NA, N,M, V,W) 

CALL MMUL(NA,m,IA,N,N,M,U,W,X) 

CA LL TRNAIB (NA , M, N, M, U, W) 

CALL MMUL (NA,1SR, M, N,N,M, V,W, Y) 

CALL MADD (NA, NV, N, N, X, Y, X) 

CALL MSCAIE (NA,N,N,-1. ©0, X) 

CALL TRNATA (NA, N, X) 

CALL WEIGHT NAA, ICON, NA, N, N, P (1, K) ,S, V) 

CALL MAED pIA, NA, tA , N, N, X, V, X) 

CALL WEIGHT (NB,NAA, ICON, KA,M,M, P (1, K) ,BS8, Y) 

CALL MAEO(NA,NA, Ni,M,M,Y,R, Y) 

CALL MMUL (NA,bA, NA,M, N,M, U,Y,V) 

CALL MMUL(NA,NA, 1SA,N,N,M,V,W,Y) 

CALL MAED (NA, m, NA, N, N, X, Y, X) 

CALL MCF (NA, NA, NA, N, N, X, A, U, WORK) 

CALL MAED (NQ, NA, NS , N, N, Q, U, S (1, 1, K) ) 

1000 CCNTINUE 

IF (ICOUNT-IEND) 11,12,12 

11 ICOONT= ICOUNT + 1 
IF(ICOUNT.IT.IPRT) GO TO 1 
ICIMl = ICOUIB? -1 

WRITE (ROOT, 5000) ICIMl 

DO 1005 k=i,k:cn 
ml = K-1 

WRITE (KOUT, 4000) KMl 

CALL MATIO (NS, N, N, S (1 , 1, K) , 3) 

1005 CONTINUE 
GO TO 1 

12 CONTINUE 
C 

C CCMPUTE OPTIMAL COST FENCTION 

CAIL WEIGHT (NA,NAA, ICON, NA,N,N,E,S,U) 

WRITE (ROOT, 7000) 

CALL MATIO (NA, N, N, U, 3) 
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GOTO (23, 22),JEND 


CCMFUTE G OPT 

23 CALL MM^JL{NA,^©l,m,N,M,N,W,A, U) 

CALL M9CALE (NA,N,N,-1. 0D0,U) 

WRITE TKOUr, 6003 ) 

CALL MATIO (NA,M, N, 0, 3) 

IX) 217 KP=1,BCCN 

CALL MMUL(NA,^B,^A,N,N,M,B(1, 1,KP) ,U,W) 
CALL MAED(NA,NA,ffi,N,N,A,W,W) 

CAIL EIGVAL(NA,N,W,W,FR,FZ,WGRK,IPVT) 
,217 CCOTINOE 
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CAICULATE CCMIARISaSI WITH GNORM 
ICOUNT= 0 
DO 130 K=1,ICCN 
DO 120 J=1,N 
DO 110 1=1, N 

S (I, J,K) = 0. 0D0 
110 carriNUE 
120 S(J,J,K) = l.©0 
130 CONTINUE 
JEND= 2 
400 CONTINUE 

DO 98 K=l,K:a!J 

CALL WEIGHT (NS, NAA, HOCN, ISA, N, N, P (1 , K) ,S, U) 

CALL MQF’(NA,NA,IA,N,N,U,A,X,WCRK) 

DO 96 L=1,ICC»I 

CALL Mff(NS,I©,IB,N,M,S(l,l,L) ,B(1,1,L) ,SB(1,1,L) ,WORK) 
96 Ca>mNUE 

CAIL WEIGHT (NS,MA, ICCN,m,M,M,P(l, K) ,SB, Y) 

CALL MCF (NA,NA, IA,M, N, Y, GNORM (1, 1, K) ,U, WORK) 

CALL MAID (NA, NA, KA, N, N, U, X, X) 

DO 95 l=i,r:cn 

CALL MMUL(NS,m,lS,M,N,N,S(l,l,L),B(l,l,L),SB(l,l,L)) 

95 ca>?riNUE 

CALL WEIGHT (NS, ^A, £CCN,bA, N,M, P (1, K) ,SB,Y) 

CALL TRNATB (NA, KA, N,M, Y,W) 

CALL TRNATA(NA,N, A) 

CAIL MMUL(NA,m,m,M,N,N,A, Y,V) 

CAIL MMUL{NA,NG,NA,N,N,M,V,GN0BM(1,1,K) ,Y) 

CAIL MAQ3 fNA, m, N, N, Y, X, X) 

CALL TRNATB (NG, tSA , M, N, ®ORM ( 1 , 1, K) , V) 

CALL MMUL(NA,Mi,NA,N,N,M,V,W,U) 

CAIL TRNATA (NA, N, A ) 

CALL MMUL(NA,l«^,m,N,N,N,U,A,W) 

CAIL MACD(NA,l^,m,N,N,W,X,X) 

CALL MADD (NA, m, m, N,N, X, Q, X) 

CAIL MCF (NR, NG, 1SA,M, N, R, OsIORTI (1, 1, K) ,U,W®K) 
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CALL MAE© (NA, NA, N,N, X, U, X) 

CALL SAVE (NA,NS,N,N, X,S{1,1,K)) 

98 CONriNUE 

IF{ICOUNT-IEeB) 4010,4011, 4011 
4010' ICOUNT= iCOUNT + 1 
00 TO 400 

4011 WRITE (KOUr, 9000) 

CALL MCF(NA,NA, IA,N,N, X,A,U,WCRK) 

DO 1306 l=i,k:cn 
mi = L-1 

WRIIE (ROOT, 4000 ) LMl 

CALL MATIO , N, W, S (1 , 1, L) , 3) 

1006 CCNIINOE 
<30 TO 12 

4000 FORMAT (/, 41 S,I5,/) 

5000 FORMAT {//,1IH ITERATION ,13) 

6000 FORMAT (//,10H G OPTIMAL ) 

7000 FOiMAT(//,39H OPTIMAL COST FIKCTION X'CX, WHERE C IS,/) 
8000 FORMAT (//, 31 P,/) 

9000 FOK4AT(//,38H COST CCMIARISCM WITH NORMAL SOLUTION 
9500 FOE^AT{2D25.15) 

9600 FORMAT (/,3H A ) 

9700 FORMAT (/,3H Q ) 

9800 FORMAT (/, 31 R ) 

9930 FORMAT (/, 31 B, 15,/) 

44442 FORMAT!/, 31 S ,15,/) 

44443 FCRMAT{/,13H A + B*GZERO) 

44441 FORMAT (/,45H SOLUTION TO STANDARD OPTIMAL CONTROL PROBLEM) 
2 STOP 
22 RETURN 
END 
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SUBROUTINE SWITCH {SIA,RB, bC, N3, NAR, RAC, N, IR, ISAA, ICCN,M, A, B, P,C, G, 
1X0, E, BTEMP, EM, WORK, Y, U, V,W, W, IPVT, ARRAY, DT, NPOINT, NGRIDH,MCa<l) 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*****PARAMETERS: 

INTEOIR NA,bB,lSC,bAR,lRC,N, IR,bAA,BCCN,M, NroiNT,N3 
INTE(ER r4Ca^ ^POINT) ,IPVT (N) 

DOUBLE PRECISION A(NA,N) ,B{NB,NAA, KZCN) ,C(NC,N) ,X0 (N) 

DOUBLE PRECISION G (NG,NAA, ICON) ,Y{N) , WORK (N) ,EM(NA,N) 

DOUBIE PRECISION U (M) ,W(NA, PCCN) ,W(NA,N) , VPIA,N) 

DOUBLE PEIECISION ARRAY (NAR, bAC) ,P 0;«A, ICON) ,E (KCCN) ,ErEMP{KCOsI) 

*****LOCAL VARIABLES: 

INTEGER IN (27) ,NSYvi(l},rr(10,l) 

DOUBLE PRECISION WT{10) ,SUyi,T/JOR[, YVIIN, Y4AX, SF(10) ,ZERO, X4AX,T, DT 
DOUBIE PRECISION DD 
DIMEJISION R (30) 

*****stBROUTINES CALLED: 

MMUL,MSCAIE ,MEXP, SAVE,,FIG, IHFLT 


***** FUNCTIONS; 
GGUB,UCALC 


*****PURPOSE; 

IHIS DOUBIE IRECISION SUBROUTINE PERFCEMS THE CCMIUPATIONS 
Ab® PRINTS THE DATA FCR SIMULATION OF TOE SWITCHING GAIN 
mOBLEM RELATING TO THE PUBLICATION: 'ON THE RELATIONSHIP 
BEIWEEl'j RELIABILITY AND LINEAR QUAERATIC OPTIMAL CONTROL ' 

BY J. DOUGLAS BIRDWELL AND M. ATHAbB. 

*****PAiy^MErER DESCRIPTION; 

NA, bB,NC,bC, THE FIRST DIMENSION OF TOE ARRAS A (AND EM, 

bAR W,W,V) ,B,C,G,AND ARRAY RESPECTIVELY AS 

CECIARED IN THE CAILING PROGRAM DIMENSION 
STATEMENT ; 

bAC COLUMN DIMENSION OF TOE ARRAY CONTAINING ARRAY 

AS DECIARED IN THE CAILING PROGRAM DIMENSION 
STATEMENT ? 

N NUMBER OF STATES ; 

IR NUMBER OF OUTPUTS ; 

mh THE SECQ'ID DIMQBION OF TOE ARRAS B AND G AS 

TEC EARED IN THE CAILING PROGRAM DIMENSION 
STATEMEbJT ; 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

K, 

c 

c 

c 

c 

c 

c 

c 

c 

I-t 

V- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ICCN 

M 

A 

B 

c 

G 

M 

MCCN 

£ 

ETEMP 

WCRK 

Y 

U 

V,W,VW,EM 

IPVT 

ARRAY 

DT 

NPOINT 

NGRIDH 


THIRD DIMENSION OF THE AREA'S B AND G AS 
EECIARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT ; 

NUMBER OF CavITROLS ; 

N BY N SSTBM MATRIX? 

N BYM BY KCCN SET OP OUTRJT MATRICES; 

IR BY N OUPFUT MATRIX; 

M BY N BY KCCN SET OF FEEDBACK MATRICES.; 

INITIAL CONDITION VECTOR OF lENCJTH N ; 

VECTOR OF LENGTH NPOINT CONTAINING THE.EJACT 
COTFIGURATION INDICES,; 

XRATCH VECTOR OF LENGTH KCCN; 

SCRATCH VECTOR OF lENGTH KCCN ; 

SCRATCH VECTOR OF lENGlH N ; 

VECTOR OF lENGTH N ; 

VECTOR OF lENGTH M ; 

N BY N SCRATCH ARRAS ; 

SCRATCH VECTOR OF LENGTH N; 

NAR BY NAC WORKING ARRAY; 

NAR MLBT BE GREATER THAN OR EQUAL TO NSTEES + 1 
NAC MIBT BE GREATER THAN OR EQUAL TO IR + M ; 

STEP SIZE; 

NUMBER OF STEPS +1; 

NmBER OF I«IAJQR ORDINATE DIVISIONS USED 
IN HOTTING 

NGRIDH MLBT BE lESS THAN OR EQUAL TO 12; 


*****NOTES: 

BOTH THE OOTHJT AND THE CONTROL U(T) = -G (I)*X fT) ARE CCMRJTED, 
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C 

C GGUB IS A RANOa^ NUMBER GENERATOR 

C 

C CEAIC IS A LBER-SUPPLIED, APPLICATION SIECIFIC FCNCTION TO 
C CAICUATE THE Ca'iTROL U. 

C 

C *****HISTORY: 

C WRIFTEN BY J. A. K. CARRIG {ELEC. SK. lAB., 35-307, 

C CAMBRIDGE, MA 02139, ffl. : (617) - 253-2165), JANUARY 1978. 

C MC6T RECEMr VERSION: MARCH 22, 1978. 

C 

C 

C 

COMMCNANOUAIN, bout 
LCCN = 1 

DATA YSF/10*1.©0/,IBIANK/4H / 

DATA TWOPI/3. 1459/ 

DATA M9C,MA)ES, IXY, lEGY, ZERO,MM,NEG, IZERO/1, 0, 0, 1, 1. ®0, 1, 0, 0/ 
DATA IK(1),IN(2),IN(3),IN(4)/4H1 , 4J2 ,4H3 ,4H4 / 

DATA IN (5) ,IN{6) ,IN(7) ,IN(8)/4H5 ,4H6 ,-Si7 ,418 / 

DATA IN{9) ,IN(10},IN{ll),IN(12)/4H9 ,4110 ,4111 ,4112 / 

DATA IN{13),IN(14),LN(15),IN(16)/4H13 ,4H14 ,4H15 ,4116 / 

DATA IN(17),IN{18),INa9),IN{20)/4H17 ,4118 ,4119 ,4120 / 

DATA IN (21), IN (22), IN (23) , IN (24) /4H 21 ,4122 ,4123 ,4124 / 

DATA IN (25), IN (26), IN (27) /4H 25 ,41 Y, 41 U/ 

DATA IT (3, 1) ,rr (4, 1) ,11 (5, 1)/4HVER3, 4UB T, 41IME / 

DATA 11(6,1), IT (7,1), 11(3, D/4H ,41 ,41 /' 

DATA IP(9,1),IT(10,1)/4H , 4H / 

IX=35 

DO 61 IZ=1,NP0INT 
61 MCCN (IZ)= MCON(IZ) + 1 
‘IWCH = 2. £D0*rWOPI 
NSTEPS ■= NPOINT -1 
T= 0.fD0 

3031 FORMAT (24H EJACT CCNFIGURATION = ,13) 

CALL MMUL (NC , N, N, MM , IR, N, C, X0, Y) 

CALL MMUL(NA,N,M,MM,M,N,G(1,1,MCCN (1)} ,)^,U) 

WRITE (KOUr, 1500) 

WRllE (KOUr, 1200) 

WRITE (KOUT, 1300 ) 

WRITE (KOUr, 1000) T 
1031 FORMAT (/,12H GAIN MATRIX) 

WRITE (KOUT, 1100 ) (Y (I ) , 1=1 , IR) 

WRITE rxour, 1102 ) (u (i ) , 1 = 1 , m) 

C WRITE (KOUT, 1001 ) 

30 ARRAY{1,J)= Y(J) 

DO 40 J=1,M 
40 ARRAY(1,IRHJ)= U (J) 

50 DO 100 K=1,NSTEPS 
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WRITE (KOOT, 1002) K 
IF(N.EQ. 1) GO TO 72 
CALL GGLB(IX,1,R) 

WT‘(2)= TWOPI*R-(L) 

CALL GGUB(IX, 1,R) 

WT(1) = R(1)*DCCS (WT(2)) 

Wl (2)= R(1)*DSIN (WT(2)) 

GO TO 73 

72 CALL GGU3(IX, 1,R) 

WT(1) = (R (1)*2. ®0)-l.ffi0 

73. CALL MI*1UL(NA,N,N,MM,N,N,EM,WT,W®K) 

CALL MMUTi(NA,N^N,wiM,N,M,B(l, 1,MC0N (K)) ,U,BrBMP) 
CALL MAEOC!J,N,N,N,MM,ErEMP,WCRK,ETEMP) 

CALL MMUL (NA,N, N,MM, N, N, A, X0,WORK) 

CALL MADD {.'I, N, N, N,MM, ETE!4P,WCRK, X0) 

DO 52 KK = 1,ICCN 

CALL MMUL(NA,N,N,MM,N,M,B(1,1,KK),D,Y) 

CALL MSlB(N,N,N,N,MM,ErEMP,Y,Y)' 

SIM= 0.0D0 
DO 55 IIJ= l,ICON 

55 SIM = Sm + Y(IIJ)*y (IIJ) 
sm = DSCRT (SIM) 

Wl'(KK) = 0.2D0 

56 IF {StM-I£. 1.000) WT(KK) = 1. 0D0’ 

52 CONTINUE 

CALL FIG (KCCtSi,E,£TEMP,WT,ICQ'I) 

831- FamT(18H PEfT-l/T-l) = ,€>25.15) 

WRITE (ROOT, 881) (E (10) ,10=1, ICON) 

CALL MMUL (NA, ICON, ICON,!, ICON, ICON,?, ETEMP,E) 
WRITE (KOUT, 882) (ETEMP(IO) ,10=1, ICON) 

882 FORMAT (18H PI (T-l/T) = ,4)25.15) 

LCCNMl = LCCN -1 
WRITE (KOUT, 4001) LCONMl 
MCCNMl = MCON (K+1 ) - 1 
WRITE (KOUT, 3001 ) MCONMl ' 

4001 FORMAT (29H CAICULATED CONFIGURATION = ,I3) 
1002 FCH4AT(/,10H TIME STEP, 13) 

CALL MMUL(NC,N,N,MM,IR,N,C,X0, Y) 

CALL MI4UL(NA,N,M,MM,M,N,G(1, 1,IC0N) ,X0,U) 

DO 70 11= 1,M 

70 U(M) = UCArC(U,EM,B(l,l,l),B(l,l, 2)) 

T= T+ DT 

WRITE (KOUT, 1100 ) (Y (I ) , 1=1 , IR) 

WRITE (KOUT, 1102) (U (I) ,I=1,M) 

DO 80 J=1,IR 

80 ARRAY-{14K, J)= Y(J) 

DO 90 J=1,M 

90 ARRAY (14K, IR-0 )= U (J) 

100 CCNTINUE 
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»1AX = DFIDATiNSTEfB)*^ 

IW= ROUT 
NS^(1)= 25 
IT(lfl)= IN (26) 

DO 110 J=1,IR 

IF(J.I£.25) IT (2,1)= IN(J) 

IF (J. or. 25) IT (2,1)= IBIANK 

110 CALL THPLT (IW, lEG^, NK)INT,ZERO, )MAX,NGRIDH, MIN, MAX, MF, IT 
.1 ARRAY (1, J) ,NAR,NLG,MSC, MARES, m,NS^) 

IT (1,1)= IN (27) 

NSM(l) = 21 
DO 120 J=1,M 

IP (J. LB. 25) IT (2,1)= IN(J) 

IF (J. or. 25) IT (2,1)= IBIANK 

120 CALL THHiT (IW,IEGy,NR)INT,ZERO,X4AX,NGRIDH, Y^lIN, Y'1AX,’SF, IT 
1 ARRAYd, J+IR) ,mR,NLG,MSC, MAXES, IXZ, NS ^) 

1100 FCmAT (4H Y = , 5 (2X, 1PD19. 8) ) 

1000 FCM1AT (5H T = , F5. 2) 

1102 FORMAT (4H U = , 5 (2X, 1PD19. 8) ) 

1200 FORMAT (IIH OUTRJT Y) 

1300 FORMAT (12H CCNTROL U) 

1400 FCH4AT(/,28H SIMUEATIOSI OF LINEAR SMTEM,/) 

1500 FORMAT (/,3IH SIMULATION OP LINEAR REGULATOR,/) 

RETURN 

END 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBKJOTINE READY (NAA,NA, bB, NQ, NR, N3, NS, m, N,M, ICON, A, B, R, Q, P, 
MR,WI,S,SB, 0, V,W,X, Y,<aiORM, R^D, RADINV,BS8,HCRK, IPVT, IEND,ieTEK) 

*****PARAMETERS_: 

INTEGER NAA , m , NQ, NR, N3, NS , MIA, N, M, ICON, IPVT G'l ) 

DOIBIE PRECISION A(NA,N) ,X(NA,N) ,Q(NQ, N) ,R{NR,M) 

DOUBIE PRECISION S (NS,NAA, RCCN) ,P (NA, ICON) ,EB(NS,JSRA, ICON) 

DOUBLE PRECISION GN0RM(NG,I;AA,K:CN) ,BSB(NB, NAA, ICON) ,WR(N) ,WI (N) 
DOUBIE PRECISION B(N3,I^A, ICON) , RAD (NRA,N) ,RADINV(NRA,N) 

DOUBIE PRECISION U {NA, W) , V(NA,N) ,W(NA,N) , Y (NA,N) ,WORK (N) 

*****LOCAL mRIABLES: 

DOUBLE PRECISION CCND 

INTEGER KIN, BOUT, KL, Hyl 1, J, I, K, tIEND, NEND, L, IM 1 
*****SlBROUTINES CALLED: 

MQF ,MADD, MUNEQ, 'mNATB,MMUL,MSCAIE , EIGVAL, SAVE, WEIGHT 


******PUrPOSE: 

■IHIS DOUBIE PRECISION SUBROUTINE SOLVES THE SWITCHING-GAIN PROBIEM 
RELATING TO TOE PUBLICATION: 'ON THE RELATIONSHIP BEB^EEl^I 
RELIABLILITY AND LINEAR QUAERATIC OPTIMAL COOTROL ' 

BY J. DOUGLAS BIRCWELL AND M. ATHAIB. 

*****PARAMETER DESCRIPTION: 

ON INPUT: 

ISAA IHE SECOND DIMENSION OF THE ARRAS S,®,QSKIRM, 

BSB,B AS EECIARED IN THE CALLING PROGRAM 
DIMENSION STATEMEJJT; 

THE FIRST DIMENSION OP THE ARRAS 
A (AND P,X,U,V,W,Y) ,B{AND BSB) ,Q,R,OaORM, 

S (AND SB) ,RAD (AND RADINV) RESPECTIVELY 
AS EEC EARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT ; 

THE NUMBER OF STATES ; 

THE NUMBER OP CBSERmTIONS ; 

THE NUMBER OF CCNFIGURATIONS ; 

N BY N SSTEM mTRIX; 

N BY M BY KCa^ SET OF INPUT MATRICES; 

M BY M CONTROL WEIGHTING MATRIX; 


C«i,NB,N0,NR, 

NG,NS,ERA 

N 

M 

ICON 

A 

B 

R 
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C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

V-. 

c 

c 

c 

c 

c 

c 

/~l 

c 

c 

c 

c 

c 

c 

r> 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Q 

p 

asi OOTRfT: 

WR,WI 

s 

SB,B,BS8 
U,V,W,X, Y 
GNORH 


RAD, RADINV 
WORK 
IPVT 
lEND 


N BY N STATE WEIGHTING MATRIX; 
KCQ'>i BY KCCN PROBABILITY MATRIX; 

SCRATCH VECTORS OF lENGIH N; 

N BY N BY KCCN SET OF SOLUTIONS ; 
N BY N BY KCCN SCRATCH ARRAYS ; 

N BY N SCRATCH ARRAB ; 


ORIGINAL PAGE IS 
OE POOR QUALITY 


N BY M BY KCCN ARRAY USED TO STORE THE 
QUN MATRICES FCR THE NORMAL LINEAR (3UAERATIC 
(AIBSIAN PROBLEM. CN RETURN, GNORM CCKTAINS THE 
G^INS ASSOCIATED WITH THE SWITCHING GAIN PRCBIEM; 

N BY N SCRATCH ARRA^ ; 

SCRATCH VECTOR OF LENGTH N; 

SCRATCH VECTOR OF .lENGTH N ; 

NUMBER OF TTERATIONS USED IN SOLVING THE NORMAL 
LINEAR QUAERATIC (AIBSIAN PROBLEM; 


^BTEES 


NUI4BER OF TIME STEPS USED IN CCMFOTING S 


*****NOTES; 

THE SOLUTIONS TO THE NORMAL LINEAR C3UAERATIC HICBIEM, 

■IHE EIGEN\ALUES OP ‘IHE MATRICES (A + 3(I)*GN0Rf^(ER0) ) 

AS WELL AS THE EIGENVALUES OF THE MATRICES (A + B(I)*G(I)) 
ARE PRINTED. 


*****HISTORY: 

WRITTEN BY J. A. K. CARRIG (ELEC. SYS. lAB. , M. I.T. , EM. 35-307, 
CAMBRIDGE, MA 02139, PH.: (617) - 253-2165), JANUARY 1978. 
MCST RECENT VERSION: MARCH 22, 1978. 


ccmmcn/inouain, lout 

WRITE (KOUT, 96 00) 

CALL MATIO CIA, N, N, A, 3) 
WRITE fKOUT, 9700) 
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CALL MATIO (NA, N, N, Q, 3) 

WRITE <KOUr, 9800) 

CALL MATIO (InIR, N, N, R, 3) 

DO 222 KL=1,H:CK 
Ml = KL-1 

WRITE (ROUT, 9900) KMl 

CALL MATIO (N3,N,M,B{1, 1, KL) ,3) 

DO 4 J=1,N 
DO 3 1=1, N 

3 Y(I,J)= 0»ID0 

4 Y(J,J)= l.eD0 
DO 210 K=1,IEND 

CALL Mff N,M, Y,B(1, 1, KL) ,U,WORK) 

CALL MAIX)(NA,NR,IA,M,M,U,R, U) 

DO 14 J=1,M 
DO 13 1=1, M 

13 V(I,J)= 0.3)0 

14 V(J,J)= 1.0)0 

CALL MLINEQ(NA,m,H,M, U, V,COND, IPVT,WCRK) 

CALL TRNATB(NB,m,N,M,B(l, 1,KL),X) 

CALL MMUL(NA,m,m,N,M,N, X,Y,U) 

CALL MI4UL (NA, m, lA , N,M, N, U, A, X) 

CALL MCP (NA,IA,IA,M,N, V, X,W,WORK) 

CALL MSCAIE <NA,N,N,-1. 0)0, W) 

CALL MS’{NA,lA,m,N,N, Y,A, U,WGRK) 

CALL -MAED (NA, lA, lA, N, N, U, Q, U) 

CALL MACD fNA, m, lA, N, N, U,W, Y) 

210 Ca^TINUE 

WRITE {ROOT, 44441) 

WRITE (KOUT, 44442) 

CALL MATIO m,N,N, Y, 3) 

CALL MT4UL (NA, m, N3, N, M, M , V, X, O^ORM ( 1 , 1, KL) ) 

CALL MSCAIE (NA,M, N, -1. 0)0, GNORM(l, 1, KL) ) 

WRITE (KOUT, 6000 ) 

CALL MATIO (14G, M, N, (UORM ( 1 , 1, KL ) , 3 ) 

CALL MMUL(N3,NG,tA,N,N,M,B{l,l,KL) ,Q^ORM(l, 1, 1) ,V) 
CALL MADD (NA, KA, ^A , N, N, V, A, V) 

WRITE fROUI, 7008) 

CALL EIG\AL fNA, N, V, V,WR,WI,WCRK, IPVT ) 

222 CCSSniNUE 
JEND= 1 
26 CCWIINUE 

WRITE fiCOOT, 8000) 

CALL MATIO (NA, BCCN, ICCN, P, 3) 

DO 5 k=i,k:cn 

CALL SA\E (NQ,NS,N,N,Q,S(1,1,K)) 

5 CCNIINCJE 

DO 91 NEND= 1,NSTEPS 
WRITE (KOOT, 4500) NEND 
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20 a COWTINUE 

DO 90 L=1,£CCN 
DO 80 K=l,H:a^ 

CALL MCF (NS, bB, N,M, S (1, 1, K) ,B (1, 1, K) ,BSB( 1, 1, K) ,WCRK) 

CALL MMUL(NS,^B,^E,M,N,N,Sa, 1,K) ,B{i, 1,K),SB(1, 1,K)) 

80 CONTINUE 

CAIL WEIGHT (NS, NAA,ICa<I,m,N,M,P(l,L) ,SB,V) 

CALL WEIGHT (NS,NAA, P(1,L) ,BSB,PAD) 

CALL MAED C^1R,NRA,1SA,M,M, R,RAD,U) 

DO 98 J=1,M 
DO 97 1=1, M 

97 RADINV(I,J) = 0.0)0 

98 RADINV(J,J)= 1.0)0 INAL PAGE IS 

CALL MLINEQ(NA,NRA,M,M,U,RADINV,CCND,IPVT,WCHK) ^ POOR QUALOT 
DO 70 K=l,iCC»] 

DO 60 J=1,N 
DO 60 1=1, M 

60 B3B(I, J,K) = SB(J,I,K) 

70 CONTINUE 

CAIL WEIGHT (NS, ICCN,m,iyi, N, P(1,L) ,BSB,U) 

CALL MMUL t\IRA,Nf\, NA, N,M,M, PADINV, U,W) 

CALL MMUL (NA, m, m, N, N,M, V,W, Y) 

CALL MMUL (NA, NA, NG, N, M, N, W, A, ®10RM ( 1 , 1, L) ) 

CALL MSCAIE (NG,M,N, -1. O)0,Q^ORM(1, 1,L)) 
ml = L-1 

WRITE fl(OUr, 2005) ml 

CALL MATIO a'^G,M, N, 0JORM (1, 1, L) , 3) 

IF(NEND.NE.KBTEIS) GO TO 73 

CALL Mi'lULa'^B,NG,m,N,N,M,B(l,l,L},(lIORM(l,l,L) ,W) 

CALL MAIX)a'gA,m,IA,N,N,A,W,W} 

WRITE (KOUT, 7009 ) ml,IMl 

CALL EIGVAL(NA,N,W,W,WR,WI,WCRK,IPVr) 

73 CALL MSCAIE (NA,N,N, -1. 0)0, Y) 

CALL WEIGHT (NA,I®^, iCCN, ISA, N, N, P(l, L) ,S,W) 

CAIL MAED a'JA, ISA, lA, N, N,W, Y, Y) 

CALL M(F(NA,NA,1A,N,N, Y,A,W,WCRK) 

CALL MADD(NA,m,IB,N,N,W,Q,S(l, 1,L}) 

WRITE (KOUT, 4000 ) mi 

CALL MATIO (NS, N, N, S (1, 1, L) , 3) 

90 COSITINUE 

91 CQ'JTINUE 

2000 FCBMAT (3D25, 15) 

4005 FORMAT (3H S) 

2005 FORMAT (4H G,I3) 

4000 FORMAT (4H 3,13) 

4500 FORMAT (IIH TIME= T2 -,I3) 

5000 F0R4AT{11H ITERATION ,13) 

6000 FCBMAT(10H G OPTIMAL ) 

7000 FCRMAT(40H OPTIMAL COST FUNCTION X C X, WHERE C IS) 
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7008 FQK4AT(21fl A + B(I)*G3TAR (ZERO)) 

7009 FORMAT (7H A + B, 13, Si * G, 13) 

8000 FCmAT(3H P) 

9500 FORMAT .(3U25.-15) 

9700 FCKViAT{3H Q) 

9600 FOH^TpH A) 

9300 FORMAT (3H R) 

9900 FORMAT (3H B, 13) 

44441 FORf4AT(/,45H SOLUTION TO STANDARD OPTIMAL CONTROL PRCBIEM) 
2 STOP 

44442 FORMAT (3H S ) 

RETURN 

END 
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SlBROUriNE WEIGHT (NA, NAA, KOCN, NX, N,M, E, A, X) 


******PARAM ETERS ; 

INTEGER NA,bAA,K:CISl,NX,N,H 

EDUBIE FRECISION E (KCCN) ,A{NA,NAA, BCW) , X(NX,M) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*****L0CAL ^LABIfiS: 

INTECER I,J,K 
DOUBLE ERECISION SUM 

*****StBROUTINES CALLED: 

NONE 

ORIGINAL PAGE IS 

*****PURroSE; OF POOR QUALITY 

THIS SUBROUTINE CCMHJTES THE WEIGHTED SUM 

SlMt^TION E (I)*A(I,a,K); 1=1, N; J=1,M; K=l,fC0N. 

*****PARAMBrER EESCRIPTION: 

m THE FIRST DIMENSION OP THE ARRAY A AS DECIARED IN 

■IHE CALLING PROGRAM DIMENSION STATEMENT; 

THE SECOND DIMENSION OF THE ARRAY AS DECIARED IN 
THE CALLING PR0GRAI4 DIMENSION STATEMENT; 

ICO,'? THE THIRD DIMENSION OF THE ARRAY A AS DECIARED IN 

‘IHE CALLING PROGRAM DIMENSION STATEMENT; 

NX THE FIRST DIMENSION OF THE ARRAY X AS DECIARED IN 

CALLING PROGRAM DIMENSION STATEMENT; 

N THE ROW SIZE OF A; 


C M THE COLUMN SIZE OF A; 

C 

C E \®:CIDR OF LENGTH KCCN ; 

C 

C A N BY; M ARRAY 

C 

C *****HISTORY: 

C WRITTEN BY J. A. K. CARRIG (EIEC. SE. lAB., M.I.T., m. 35-307, 

C CAMBRIDGE, MA 02139, FH. : (617) - 253-2165), JANUARY 1978. 

C MGBT RECENT VERSION MARCH 22, 1978. 

C 

C 

C 

DO 10 J=1,M 
DO 10 1=1, N 
X(I,J) =0.0)0 



WEIGHT FCRTRAN 


DO 10 K=1,ICCN 

10 X{I, J) = X(I,J) + E(K)*A(I, J,K) 

RETURiS! 

END 



lEAIC FORTRAN 


F USlCTia^ UCAIC (0 , EM , B , C) 

DOUBIE PRECISION U (10, 2) ,EM(10, 2) ,B( 10, 2) ,C(10, 2) 

RETURN 

END 
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FIG FC»TRAN 

SlBROUriNE F IG (KCQ^, E, EirEMP,WCRK, ICCN ) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*****PARAMETERS; 

DDUBIE PRECISION WCRK.(KCCN)..,E (KCCN) ,ETa4P-(KCGN) 

*****LOCAL VARIABLES: 

INTEGER MI'I, LTEMP, IFIAG, KK, IP, lU 
DOUaiE PRECISION Stl-I 

*****3GBR0UriNES CALLED: 

NONE 


*****PURPOSE: 

THIS DOIBIE PRECISION SUBROUTINE IS USED IN HYPOTHESIS TESTING. 
AT EACH TIME T, ONE OF BCCN HYPOTHESES IS CHOSEN. 


H (T-i/T) = 
I 


RHO{XfT) - A*Xfr-l) - B(I)*U fT-l))*Pi CT-lA-1) 

I 


SU4 (RHO (X CT ) - A*X fT-1 ) -B ( J )*U fT -1 ) ) *PI fT-1 A“1 ) 

J 


HYPOTHESIS H (I) IS ASSUMED TO BE COIRECT IF 

PI CT/T-1) > PI fT-l/T) FOR ALL J NOT EQUAL I 
I J 

TIES ARE RESOLVED ARBITRARILY. 

RHO(X) DENOTES THE PROBABILITY DISTRIBUTION OF X. 

*****PARAMETER EESCRIPTION: 

CN INPUT: 

ICON THE NUMBER OF HYPOTHESES ,* 

E \ECTOR OF LENGTH KCCN CCNTAINING PI fT-l/T-1); 

WCRK VECTOR OF LENGTH RCCN CCNTAINING 

RHO(X{T) - A*X(T-I) - B(I)*U(T-1)); 


CN OUTHIT: 
ETpiP 

ICON 


VECTOR OF LENGTH RCCN TO STORE PI (T/l-l ) ; 
INDICATES WHICH HYPOTHESIS HAS BEEN CHOSEN ; 


*****HIS‘IORy! 



n o n o o n 


FIG FORTRAN 


244 


WRITTEN BY J. A. K. CARRIG (ELEC, SYS. lAB. , M. I.T., RA. 35-307, 
CAMBRIDGE, MA 02139, Hi.: (617) - 253-2165), JANUARY 1978. 
MCBT RECEOT VERSION MARCH 22, 1978. 


CaMCN/INOU/KiN, KDUr 
MM = 1 
LTEMP = DCCN 
StM = 0. 0D0 
DO 10 IP = l,iCCN 
10 sm = sm + woRK(ip)*E (IP)- 
DO 20 IP=1,ICCN 

20 ETEcMP(IP) = WORRfIP)*E (IP)/Sm 
DO 60 KK = 1,ICQJ 
•IFIAG = 0 
DO 89 IU= 1, ICON 

IF(KK.EQ.IU) GO TO 79 

IF(ETEMP{KK).GT.E(IU)) IFIAG = IFIAG + 1 
79 COWriNUE 
89 CONTINUE 

IFIAG = IFIAG + 1 
IF (IFIAG. EQ.ICCN) LCOSi= KK 
60 Cafl'INUE 

IF (ICON. EQ. 0) ICON = LTEMP 

RETURN 

END 
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C lATEST \;ERSI0N 2/9/11 

DOUBLE JREC IS ION COND,BEE,WR(10) ,WI (10) 

DOUBIE ERBCISION A(10, 3),X(10,3) 

INTEGER (2).,HKvl^BC(2.).,VPIME.(2.).,K^IME(2-) 

DOUBIE FRBCISION OSlORM (10, 3, 4) 

DOUBIE PRECISION BSB(10,3,4) 

DOUBIE PRECISION S (10, 3, 4) ,P(10, 4) ,SB(10, 3, 4) 

DOUBIE HIECISION S3T(10, 3 ) ,Q(10, 3),R(10, 3) ,B(10, 3,4) 
DOUBIE PRECISION PR (4) ,P1, P2, 12 (4 ) ,H) (10, 4) ,ES (4) 

INTEGER IPVT(10) 

DOUBIE PRECISION AZERO,Aa>3E, ATWO 

DOUBIE PRECISION RAD (10, 3) , RADINV(10, 3) ,U (10, 3) 

DOUBIE PRECISION V(10, 3),W(10 ,3),Y(10, 3) ,SCK,WORK (10) 

Ca4MCN/INOD/KIN, KOUT 

NAA= 3 

AIWO= -3. 0D0 

AZERO = -4.£D0 

AONE = 6. ffi)0 

Pl= .05D0 

P2 = .75D0 

KIN= 5 

KOUL’= & 

N= 3 
M= 3 
N2 = 6 
KCai = 3 
NS= 10 
IPRT= 17 
IEND= 25 
ICOUNT = 0 
NSTEPS = 25 
NA= 10 

NRA= 10 
NR= 10 
NB= 10 
NQ= 10 
NG=10 

22 IF (ICOUNT. NE.0) READ (KIN, 9500, END=2 ) (PR (I ) ,FZ (I ) ,I=1,N) 
9500 FCRMAT(3D25. 15) 

DO 11 JK= 1,N 
DO 11 JL = 1,N 
Q(JL,JK) = 0. 010 
R(JK,JL) = 0.010 

11 A(JL, JK) = 0. 030 
BEE = -10. me 
P(l,l) = 1.0D0-P1 
P(2,2) = 1.0D0- P2 
P(3,3) = 1.030 
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P(l, 2) = 0. 0D0 
P a, 3) = 0. ®0 
P (2, 1} = PI 
P (2, 3) = 0. 

P(3,l) = 0.0D0 
P (3,2) = P2 
A(l,l)= 8_, £D0 
A(2, 2)= 0. £D0 
A (3, 3) = -AZERO 
A(l,2) = 1. 0D0 
A(2,3) = 1..0D0 
A (3,1) = -ATMO 
A (3, 2} = -AONE 
Q(l,l)= 3.030. 

Q(2,2) = 3.0D0 
Q(3,3) = 3.0D0 
R(l,l)= l.a>0 
R(2,2) = 1. 030 
R(3,3) = 1.030 

1)= 00. 030 
B(2,2, 1) = 0.0D0 
B(2,l,l)= 9.030 
B(l,2, 1)= 0.030 
B(l, 3, 1)= 0.030 
B(2,3, 1} =0.030 
B(3,3, 1) = 1. 0D0 
3(3, 1,1) = 1.030 
B(3,2, 1) = 1. 030 
B(l,l, 2) = 0.030 
B(2,2, 2) = 0. 030 
B{2,1,2) =0.030 
B(l,2, 2) = 9. 030 
B(1,3,2) = 0. 030 
B{2,3, 2) = 0. 0)0 
B{3,3,2) = BEE 
B(3,l,2) = 1. 030 
B{3,2,2) = 1.030 
B(l,l,3) = 0. 030 
B(2, 2,3) = 0. 030 
B(2,1,3) = 0. 030 
B(l,2,3) = 0.030 
B(l,3, 3) = 0.030 
B(2,3, 3) = 0.030 
B(3,3, 3) = 0.030 
B(3,l,3) = 1.030 
B{3,2, 3) 1.030 

PR(1) = .05D0 
PR (2) = .75D0 
P(l,l) = 1.030 - PR(1) 
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P«2,l) = PR(1) 

P(3,l) = 0.0D0 
P(l,2) = 0. 0D0 
P(l,3) = 0. 0D0 
P(2,2) = 1.3D0 - PR (2) 

P(3,2)= PR (2) 

P (2,3) = 0. (D0 
P(3,3) = 1.030 

C CAIL T^ME(^mMyR,HR^lN,^C,V^mE,R^IME) 

CALL READ2 (NAA, ^A, KB, NQ, NR, N3, NS, NRA, N,M, ICCN, A,B, R, 0, P, 

1 WR,WI, S, ^ , U, V,W, X, Y, GNOH'!, R\D, mDINV,Bffl,WCRK, IPVT, lEND, 

2 NSTEES) 

2 STOP 

END 
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C lATEST VERSION 2/17/78 

DOCBIE fSEClSION E (4) ,ETEMP(4 ) ,SM, SIGMA, SIQMl, ESINV, ESIGM.A, Sn4W 1 
DOUBIfi PE^ECISION COND,UJDOLE,ri®lD3V,DOLFMl,EM(10, 2),X0(10),DIN\M1 
DOUBLE FRBCISION ARRAY(100, 50) , Y0 (10) ,U0 (10) 

DOUBLE DECISION A(10, 3),C(10, 3) ,HU, m2, X(10, 3) 

E (1) = 1.270 

EM(1,1) = 1.0D0 
E;4(2,2) = 1.070 
E (2) = 0. 0D0 
ETEMP (2 ) = E (2 ) 

E (3) = 0 . 070 ORIGINAL PAGE IS 

ETE(4P(3) = E (3) OF POOR QUALirY 

ETEMP (1) = E (1) 

E (4) = 0. 070 

DOUBLE PRECISION GNORM (10, 2, 4) 

DOUBLE HIBCISION BSB{10,4, 3) 

DOUBLE H^IBION S (10, 3, 4) ,DT, P (10, 4) ,SB (10, 3, 4) 

DOUBLE PRECISION WR (4) ,WI (4 ) ,HH (4, 4) ,XX (4, 4) ,ACL(10, 3) 

DOUBLE PRECISION SBT (10, 3 ),Q(10, 3) ,R(10, 3),B{10, 2, 4) 

DOUBLE PRECISION PR (4) ,,FZ(4) ,PD (10, 4) ,P5 (4) 

INTEGER IPVT (10) ,MCCN (100) ,NSTEFS, NSRIDH, ICON (100) 

DOtJBIE PRECISION RAD(10, 3) ,I?VDINV(10, 3) ,SNH«7(10, 3, 4) ,U (10, 3) 

DOUBLE PRECISION V(10, 3) , W(10, 3) ,W(10 ,3),Y(10, 3) ,StM,WCRK (10) 

LOGICAL NOISE 

CCMI4CN/TNOUAIN, K7Ur 

K0=0 

lA = 1 

READ (5, 11111) NPOINT 
33333 READ (5, 11111, END=22222) ITIME,K 
Hill FCRMAT(2I4) 

DO 44444 IXYZ = lA, ITIME 
44444 MCCN (IXYZ) = K0 
MCCN (ITME)= K 
lA = ITIME 
K0= K 

GO TO 33333 

22222 DO 55555 IXYZ = ITIME, NPOINT 
55555 MCQ4 (im) = K0 

LUDOLP= 2. 718281828459045D0 

LUDINV= 1. 070/LUDOLF 

DOLFMl = LUDOLF - 1. 070 

DINVMl = LUDINV - 1. 070 

NAA = 2 

NC = 10 

KIN= 5 

KOOT= 5 

N= 2 

M= 2 

N2 = 4 
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k:cn = 4 
NH= 4 
N3= 10 
IPRT= 17 
IEND= 50 
IPRT =49 
icouNT = a 
NA= 10 
•W3= 10 
N>1=^ 

NRA= 10 
NR' 10 
NQ= 10 
NG=10 

PZ(1) = ,]D0 
PR{1)= .B0 
DO 15 1=2, N 
PR(I)= PR(1) 

15 PZ(I)= PZ(1) 

22 IF(ICOU'IT.NE.0) READ (KIN, 9500, END=2) (PR (I ) ,FZ (I ) ,I=1,N) 
SI®IA= 1, (D0 
ESIGMA= LUD0LF**SIQ4A 
ESINV= LUDINV**SIGMA 
C(l,l) = 1. 0D0 
C(2, 2) = 1. 0D0 
C(l, 2) = 0. £D0 
C(2,l) = 0. ©0 
or = 1. ©0 

NSTEES = 50 
A(l,l)= ESIGMA 
NAR= 100 
NAC = 50 
A (2,2)= ESINV 
A(2,l) = 0. ©0 
A(l,2) =0.©0 
Q(l,l)= 14. ©0 
Q(2,l)= 8.©0 
Q(l,2) = 8. ©0 
Q(2, 2) = 6. ©0 
R(l, 1)= 1. ©0 
R(2,l) = 0. ©0 
R(l,2) = 0. ©0 
R(2, 2 ) = 1. ©0 
B (1,1,1)= ESIcm -l.©0 
B(2,l,l)= ESINV- l.®0 
B(2,2,l) = -B(2,l,l) 

B(1,2,1)= B(1,1,1) 

B(l,l, 2) = 0. ©0 
B(2, 2, 2) = -DINV41 
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B(2,l, 2)= 0.0D0 
B(l,2, 2) = DOLFMl 
8(1,1, 3) = DOLFMl 
B{1, 2, 3) = 0. ©0 
B(2, 1, 3) = DIN\iMl 
B(2,2, 3) = 0. ©0 
PRl = .1D0 
PR2= . 1D0 
P{1,1) = .81D0 
P(2,2) = .09D0 
P(3,2) = 0. 09D0 
P(3,3) = .09D0 
P(l,2) = .81D0 
P(3,l) = ,09D0 
P(2,l) = .09D0 
P(l,3) = .81D0 
P(2,3) = 0.09D0 
P(4,l) = .0100 
P(l,4) = .81D0 
P.(4,2) = .01D0 
P{4,3) = .01D0 
P(4,4) = ,01D0 
P(2,4) = .09D0 
P(3,4) = .09D0 
WRITE (ROOT, 9903) 

CALL MATIO fi^JA, PCCN, H3CN, P, 3) 

C WRIfE (KO[Jr,46) 

46 FORMAT (/, 41 PI,/) 

47 FCK4AT (3D25. 15) 

WRITE (ROUT, 9600) 

CALL MATIO ?S1A,N, N, A, 3) 

WRITE fKOOr, 9700) 

CALL MATIO gSIA,N, N, Q, 3) 

WRIIE (ROUT, 9800 ) 

CALL MATIO a^R, N, N, R, 3) 

DO 222 K=1,BCCN 

K41 = R-1 

WRITE (KOUr, 9900) KMl 

222 WRITE (ROUT, 9500) ( (B(I, J, K) , J=1,M) ,I=1,N) 

44 CONTINUE 

DO 14 IN=i,50 
DCON (IN) = LCCN (1) 

14 CCNTINUE 

667 FORMAT (515) 

X0(1) = .02D0 

GN0RM(1,1, 1)= -1. 063361 84D0 
GNORM (2, 1, 1)= -7. 9015188D-1 
GN0R!4(1, 2, 1)= -1. 88787889D-02 
GN0RM(2, 2,1)= -5.83582495D-02 


0?pm'^ page is 
Poor quality 
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r* 

C 

c 

c 

c 

c 

c 

c 


GNORiM(l,l, 2)= 
QJ0Ri4(2, 1, 2)= 
GN0RI4{1, 2, 2)= 
(SI0R!«1(2, 2, 2)= 
GNORf^(l,l, 3)= 
Gt«0RM(2, 1, 3)= 
GN0RM(1,2, 3)= 
GN0RM(2, 2,3)= 


-3. 69012096D-31 
-1. 14016534D0 
1, 04948339D-01 
-1 . 363087-670-31 
-1. 4256676700 
-2. 874513080-31 
1. 518842850-32 
-7 . 270124380^32 


IR = 2 

NPRFL = 1 

DO 57 IK = l/iCavi 


IKMI = IK - 1 


WRITE (KOUT, 9902) iml 

57 WRITE (KOUT, 9500 ) ( (GNOR'4 (I J, IL, IK) , IL=1 , N) , I J=1 , N) 
NGRIDH = 5 
Va,l) = B(l,l,l) 

V(2,2) = B(2,2,l) 

V{2,1) = B{2,1,1) 

V(i,2) = B(l,2,l) 

CALL MI4UL (NA, NA, m, N, N,M, V, Q40RM, U) 

CAIL MAK) e^A, N, N, U, A, ACL) 

CALL MSCAIE (NG, N,M, -1. 000, QIOFM) 
lONE = 1 

CAIL MI4UL (NC , N, N, lOTE , IR, N, C, X0, Y0) 

65 F0R4AT (IX, 30 25. 15) 

CAIL ORGSIM {NA, W3 , NG, NAR, lAC , N, IR, M, AC L, C, Q^0R4, X3,WCRK, 
lY, U, IFVT, ARRAY, DT, NSTERS , NFREL ) 

CAIL READY2 {NAA, lA, NB, NQ, NR, NG, NS, M^A, N, M, tCON, A, B, R, Q, P, 

1 WR,WI, S, SB, U, V,W, X, y, QIOBM, BAD, iADINV,BSB ,WCRK, IPVT, lEND) 
DT = 1. 0D0 


X0(1) = ,0200 

C CAIL MSCAIE (NG, N, M, -1. 0D0, OlORt-I) 

X0 (2) = 0. 0D0 

CALL SWnC H (NA, BB , N: , NG, mR, tAC , N, m, NAA , HCCN,M, A, B, P, 

1 C , QK)RM, X0, E, BTEMP, EM,WCRK, Y0, U0, V,W, W, IPVT, ARRAY, OT, NSTEFS 

2 NGRIDH, MCa^) 

9530 FCK4AT(2D25, 15) 

2000 FCS54AT(/,S>25. 15) 

9600 FORMAT (/,3H A ) 

9700 F0Rr4AT(/,2I Q ) 

9800 FORMAT (/,3H R ) 

9900 FORMAT (/,B B ,15,/) 

9933 FCRJ4AT{/,3H P ) 

9932 FC»tMAT(/,3H G ,15,/) 

2 STOP 
END 
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a^i4AT2 FC^TRZ^N 
C LATEST VERSION 2/17/78 

DOUBLE PRECISION E (4 ) ,BTEMP (4 ) ,SlK, SIGMA, SICMl, ESINV, ESIC24A, SINV^il 
DOUBLE PRECISION CCND, LUDOLF, LUDINV,,DOLFMl, EM(10, 2) ,X0 (10) ,DINW1 
DOUBLE PREClBiavl ARRAY(100, 50) ,Y0 (10) ,U0 (10) 

DOUBLE PRECISIOISI A{10, 3) ,C(10, 3) ,FRl, PR2, X(10, 3) 

E (1) = 1. I2D0 
EM(1,1) = l.gD0 
EM (2, 2) = 1.ED0 
E (2) = 0. (D0 
ETEMP(2) = E (2) 

E (3 ) ~ 0. 0- 

ETEI4P(3) = E (3) ORIGINAL PAGE Is 

ETmP(l) = E{1) OP POOP. C'TtatttV 

E (4) = 0. ®0 

DOUBLE PRECISION ®IOR1«I(10, 2, 4) 

DOUBLE PRECISION BSB(10,4, 3) 

DOUBLE PRECISION S (10, 3, 4) ,DT, P (10, 4) ,SB(10, 3, 4) 

DOUBLE PRECISION WR(4) ,WI (4 ) ,HH (4, 4) , XX (4, 4) , ACL (10, 3) 

IXIUBIE PRECISION SBT (10, 3 ),Q(10, 3),R(10, 3),B(10,2,4) 

DOUBIE I^^ISION PR (4), K (4), ro (10, 4), PS (4) 

INTEGER IPVT (10) ,MCCN (100) ,NPOINT, NGRIDH, ICON (100) 

DOUBIE PRECISION RAD(10, 3) ,RADINV(10, 3) ,SNEW(10, 3, 4) ,U (10, 3) 
DOUBIE HiECISION V(10, 3) , W(10, 3) ,W(10 ,3),Y(10, 3) , SIM, WORK (10) 
LOGICAL NOISE 

ccmmcn/inooain, lOur 

K0=0 
lA = 1 

READ(5,lllll)NPOINT 
33333 READ (5, 11111, END^22222)ITIME, K 
11111 FORMAT (214) 

DO 44444 IXYZ = lA, TTIME 
44444 MCCN (IXYZ)=K0 
MCCN (ITIME)=K 
lA = ITIME 
K0=K 

GO TO 33333 

22222 DO 55555 IXYZ = ITIME, KFOINT 
55555 MCCN (IXYZ) = K0 

LUDOLF= 2. 718281828459045D0 

LUDINV= 1. 0D0AUDOLF 

DOLFMl = LUIX)LF - 1. ®0 

DINVMl = LUDINV - 1. 03 0 

NAA = 2 

NC = 10 

KIN= 5 

KOUT= 6 

N= 1 

M= 1 

N2 = 2 
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KCCN = 2 
NH= 4 
NS= 10 
IPRT= 17 
IEND= 50 
IPRT =49 
ICOULSET = 0 
NA= 10 
N3= 10 

mA= 10 
NR= 10 
NQ= 10 
NG=10 

PZ{1) = .1D0 
PR(1)= .3D0 
DO 15 1=2, N 
PR(I)= PR(1) 

15 PZ(I)= PZ(1} 

22 IF(ICOUNT.NE. 0) R£AD (KIN, 9503, END=2) (PR (I) ,FZ(I) ,I=1,N). . 
SIGMA= 1.030 
ESIOlA= LUD0LF**3ia4A 
ESINV= LUDINV**SICMA 
C(l,l) = 1. ®0 
C(2,2) = 1.330 
C(l, 2) = 0.330 
C(2,l) = 0.0D0 
DT = 1. 0D0 
A(l,l)= 1.4140D0 
NAR= 100 
NAG = 50 
Q(l,l) = 3.030 
R(l,l)= 1.030 
R(2,l) = 0, 330 
B(l,l,l)= 2.030 

B(1,1,2)= .330 
P(l,l) = .7D0 
P{2,2) = .7D0 
P(3, 2) = 0. 09D0 
P(3,3) = ..09D0 
P(l,2) = . 330 
P (3,1) = . 09D0 
P(2,l) = . 330 
P(l, 3) = .81D0 
P(2,3) = 0. 09D0 
P(4,l) = .01D0 
P(l, 4} = .81D0 
P(4,2) = . 01D0 
P (4, 3) = . 01D0 
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C 


P{4,4) = .01D0 
P(2,4) = .09D0 
P(3,4) = .09D0 


WRITE {KOUr, 9903) 

CALL MATIO a'^A, iCCN, iCCN, P, 3) 

WRITE (KOUf, 46) 

46 FGRyiAT(/,4i PI,/) 

47 FORMAT (3D 25. 15) 

WRIIE fKOUr, %00) 

CALL MATIO (MA,N,N,A, 3) 

WRIIE (KOCr, 9700) 

CALL MATIO (NA, N, N, Q, 3) 

WRITE (KOUr, 9800 ) 

CALL MATIO (NR, N, N, R, 3) 

DO 222 K=i,s:asi 

KMl = K~1 

WRITE (KOUT, 9900) KMl 

222 WRITE (KOUT, 9500) { (B(I,a, K) , J=1,M) ,I=1,M) 
44 COMIINOE 


ORIGINAL PAGE IS 
OF POOR QUALOT 


®10RI4( 1,1,1) = -1. 063361 84D0 
GNOR4(l,2,l) = -1. 88787889D-02 
O10Rf4(2,l,l) = 7 . 90 1 51 8840-^1 
GN0RM(2,2, 1} = -5. 835824 6D -02 
GlN0RI'1(1,1, 2) = -3. 69012096D-01 
GN0H4{1, 2, 2)=1. 04948339D-01 
GsJORM(2,1,2) = -1. 140163 54D0 
GM0RM{2 , 2, 2) = -1. 36 3087 67D-01 
GN0RI4 (1,1,3)= -1. 42566767D0 
GNORI4(2,l, 3)= -2.87451308D-01 
GN0RM(2 , 2, 3) = -7 . 270124 38-02 
GNORI'i(l,2,3) = 1.518842850-^2 
GMORM (1,1,4) = 0.£D0 
GMOKM{2,2, 4) = 0.0D0 
GN0RM(1,2, 4) = 0.gD0 
GN0RM(1,2, 4) = 0. ED0 
DO 14 IN=1,50 
ICON (IN) = LCCN (1) 

14 CCMIINUE 
667 FORMAT (515) 

X0 (1 ) = . 02D0 

IR = 1 

NPRPL = 1 

DO 57 IK = 1,H:CN 

IKMl = IK - 1 

WRIIE (ROOT, 9902) liMl 

57 WRITE CKOUI, 9500 ) ( (GNORM (I J, IL, IK) , IL=1 , N) , I J=1 , N) 
NGRIDIi = 5 
V(l,l) = B(l,l,l) 

V (2,2) = B(2,2,l) 
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V(2,l) = B(2,l,l) 

V(l,2) = B(l,2,l) 

CALL MI-lUL (NA, NA, N, N,M, V, GSIOPiM, U) 

CAIL -MADD (NA,KAr NV, N, N, U, A, AC L) 

CALL MSCALE (L\1G,N,M,-1. 
lONE = 1 

C CALL MHUL (NC, N, N, lOSE, IR, N, C^ X0, !^0) 

66 FOR'IAT (IX, 3D 25. 15) 

C CAIL DR(SIM(NA,NC,NG,mR,IAC,N,IR,M,ACL,C,QsIORM,X0,WORK, 

C lY, U, IPVT, ARRAY, DC, NPOINT , NPRFL ) 

C CAIL READY2(NAA,IR,m,NQ,NR,N3,NS,NRA,HM,K:CN,A,B,R,Q,P, 

C 1 f/JR, W I, S, , U, V,W, X, Y, QiORA, lAD, mOENV, BSB ,WCRK, IPVT, lEND ) 

DT = 1. £D0 
X0(1) = -.0200 

C CAIL MSCAIE C.'^G,N,M,-1. 0D0, QJORM) 

X0(2) = 0.0)0 

CALL SWIICH(^^A,^8,N:,NS,^^\R,^AC,N,IR,KAA,H:CN,M,A,B,P, 

1 C , aORM, X0, E, BTEHP, EM,WCRK, Y0, U0, V,W, W, IPVT, ARRAY, DT, NPOINT, 

2 NGRIDH,MCQN) 

9500 FORMAT (2D 25. 15) 

2030 FORMAT (/, 3D 25. 15) 

9600 FCH'IAT(/,3H A ) 

9700 FCeHAT{/,3! Q ) 

9300 FORMAT (/, 31 R ) 

9900 F0K«1AT(/,31 B , I5,/) 

9903 FORMAT (/, 31 P ) 

9902 FCR4AT(/,31 G ,15,/) 

2 STOP 
END 
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